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The design stages of synthesis, analysis and optimization constitute the core of any 
engineering design methodology. The concept of optimiztion is intrinsically tied to hu- 
manity’s desire to excel. Though we may not consciously recognize it and though the 
optimiztion process takes different forms in different fields of endeavour, this drive to do 
better than before consumes much of our energy. Optimization is the process of maxi- 
mizing a desired quantity or minimizing an undesired one. The process of mechanical 
design, or synthesis, may often be viewed as an optimization task. In optimization, the 
goal is to minimize an objective, such as weight or peak stress, subject to constraints on 
strength and performance. Optimiztion problems may be classified into the following 
categories: 


(i) Standard or Size optimization. 

(ii) Combinatorial optimization. 

(iii) Shape optimization. 

(iv) Topology optimization. 
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Size optimization relates to determining values of design parameters such as thickness 
and cross-sectional areas and inertias. Whereas in combinatorial optimization, the prob- 
lems considered are called discrete optimization problems. In such problems some or all 
the decision variables are restricted to assume values within specified discrete sets. In 
combinatorial problems the optimum combination out of a possible set of combinations 
has to be determined. Some applications of this class of optimization include: airline 
crew scheduling, metal trimming and nesting problems. In trim problems, sheets of 
various sizes have to be cut from a big master sheet of given dimensions to minimize 
total wastage. In nesting problems, a set of component pieces in the form of specific 
shapes and sizes are to be nested on a given sheet so as to minimize the wastage of the 
resource material. Shape optimization problems deal with determininig the outline or 
shape of a component and its associated dimensions such as, for example, the shape and 
size of a fillet or the shape and size of a cutout. Shape optimization is an important 
class of structural optimization problems. The exterior and interior boundary shapes 
of structures should be controlled and determined in the shape optimization process. 
Therefore, the parameters describing the changing boundary shapes of structures will be 
taken as design variables (referred to as Shape Design Variables). Topology optimiza- 
tion addresses the basic question of designing an optimum topological structure of a 
continuum so as to satisfy a set of functional requirements and achieve the minimum or 
maximum of an objective function. Shape optimization problems deal with the design 
of the boundary while topological optimization problems deal with the nature of the 
interior continuum in terms of the distribution of its solid/void nature. 

Design optimization methods have been used successfully for optimization problems 
which involve selecting a set of optimal values for scalar design variables. The methods 
developed have been widely used for structural design as well as mechanical component 
design activities. Shape optimization is also a major consideration in many engineering 
design problems. The problem of shape optimization for a component can be described 
as deciding the geometry of a curve or a surface which will maximize or minimize an 
objective function as well as satisfy a set of constraints. The success of any shape 
optimization method depends on how the shape synthesis module has been defined. 
The most frequently published approach is that of representing the shape in terms of 
either an array of control points or a spline like representation. This representation does 
not evaluate shape explicitly and does not indicate clearly how shape should be modified 
for improving a design. Generally shape is synthesized using a mathematical model. In 
the present work, the shape synthesis is proposed to be done in the intrinsic geometry 
domain. The term shape is considered to describe the intrinsic properties of geometry 
that define a curve or surface. This representation of shape is then mapped into a 
cartesian domain, in which different analyses are performed based upon the particular 
application. The results are evaluated and updated till the optimum shape is found. 
In the last two decades, analysis using FEM/BEM (Finite Element Method/ Boundary 
Element Method) have been used for engineering applications, in general, and for the 
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analysis part of any shape optimization procedure, in particular. Numerical methods 
are needed to optimize engineering components that generally have complex loading and 
geometry. Moreover, material reduction optimization needs simulataneous consideration 
of design and manufacturing. Shape optimization has been successfully applied to a 
variety of engineering problems in the automotive, aircraft and other industries. A few 
examples are fillets, turbine disks and blades, pre-forms in forging, channels or ducts, 
minimum-stress forms for elastic solids, and air-foil shapes. 

A review of the current state-of-the-art shows that the geometrical description of a 
component is specified by using a set of auxiliary points or ai. set of points lying on the 
geometry of the curve or surface. Auxiliary points may not lie on the curve or surface. 
For example, the control polygon points of a Bezier or B-spline curve or surface are 
the auxiliary points of a particular curve or surface being designed. They, however, do 
define the shape of the curve or the surface. The review of the current literature shows 
that introducing the concept of intrinsic shape definition will definitely provide a better 
foundation for solving shape optimization problems. 

Among all the numerical techniques available for engineering analysis. Boundary El- 
ement Method is more suitable for shape optimization since one needs solutions only 
on the boundary. The shape optimization problems described above can also be formu- 
lated and solved using the shape optimization methodology based on intrinsic geometry 
proposed in this research work. The method of shape optimization, proposed here, 
consists of selecting a shape model, defining a set of intrinsic shape design variables 
and then evaluating cartesian co-ordinates of the curve. A shape model is conceived 
as a set of continuous piecewise LINear Curvature Elements (LINCEs), each element 
defined as a function of the arc length. The SDVs (Shape Design Variables) are the 
values of curvature and/or arc lengths at some of the linear curvature elements. Then, 
the mapping from the intrinsic space to the cartesian space will be carried out. Thus, 
shape synthesis will be accomplished. The shape synthesis module can be integrated 
through BEM analysis and an exhaustive search NLP (Non-Linear Programming) in 
one approach, and via first-order method of NLP in the other approach, to show how 
intrinsic geometry based shape definition can be used for engineering analysis and shape 
optimization. The proposed research also involves the method of calculating sensitivities 
of the objective functions like the frontal areas and the peak stress with respect to the 
intrinsic shape design variables. The optimal shapes of the following components have 
been obtained using the proposed shape optimization methodology: a thick cylinder 
subjected to constant internal pressure, an elastic ring under diametral compression, 
a torque arm subjected to axial and transverse loading, a fillet subjected to an axial 
load, and a ladle hook subjected to a tensile loading; The complete proposed method- 
ology has successfully been implemented as a computer aided design tool by developing 
a multi-window concept for cartesian and intrinsic space using SUN-view graphics on 
SUN-workstation with Unix environment and also using HP-9000/800 series- workstation 
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with STARBASE graphics on Unix platform. As a conclusion, it can be said that the 
proposed approach requires fewer shape design variables as compared to the methods 
where shape is represented using spline-like functions. 

The overall objective of this work is to develop an approach for shape optimization 
using intrinsic geometry and Boundary Element Method and to show how this method- 
ology can be applied to various structural and mechanical design problems. The specific 
objectives of this research work can be stated as follows: 

(1) . For the proposed methodology of shape optimization, it is required to develop 
and test a robust numerical scheme in order to solve Serret-Frenet equations. The Serret- 
Frenet equations play an important role in the curve synthesis. Since these are non-linear 
equations in intrinsic varibles, to be solved simultaneously in order to satisfy geometric 
constraints during the first stage of shape synthesis, proper caution must be taken in 
deciding free and dependent design variables. Suitable attention also must be given to 
select a stable numerical scheme and its implementation. The robust numerical scheme 
eliminates the designer intervention during the design process. This involves selection 
and development of a numerical integration module from amongst various numerical 
integration schemes such as Trapezoidal rule, Simpson’s rule and Gaussian quadrature 
etc. in order to develop a suitable scheme of shape synthesis & optimization. These 
numerical integration schemes will provide the necessary mapping from intrinsic domain 
to cartesian domain by solving the Serret-Frenet equations. 

(2) . It is also required to develop a generalized BEM analysis computer code based on 
existing sub-modules for 2D-elcistostatic plane-stress, plane-strain, and axi-symmetric 
problems which take into account the geometrical definition of the boundary in the form 
of an intrinsic curve. Boundary Element Method for 2D-elastostatic problems will be 
applied as an analysis tool in the proposed shape optimization methodology and the 
problems considered are characterized by a state of either plane-stress or plane-strain. 

(3) . The key idea is to employ geometric modeling concepts of curve synthesis typical 
of the CAD technology in order to produce sensitivity analysis results. These sensitivi- 
ties with respect to intrinsic shape design variables can then be used by an optimizer to 
generate an improved design. Presently, analysis of engineering problems requires either 
FEM/BEM approach or either algebraic or numerical analysis methods. All these tech- 
niques are carried out in cartesian domain. The intrinsic geometry should provide the 
information about the cartesian geoemetry and any change in the shape made through 
the intrinsic design variables should update the cartesian geometry. In short, it is neces- 
sary to develop a method in which intrinsic geometry along with the cartesian geometry 
will be used for engineering analysis. To study and develop the numerical aspects of 
sensitivity of the curve or boundary of the component and hence the frontal area of the 
cross-section and the magnitude of the peak stress (Maximum Von Mises stress) with 




respect to the intrinsic shape design varibles. Determination of the sensitivity of the 
design to variations of the parameters is important in design and design optimization. 
This seems particularly important since the intrinsic definition of the curve has been 
used to define the boundary of the component being designed in the present work. It is 
expected that changes in some of the intrinsic shape design variables will have a large 
impact on the resulting cartesian coordinates. 

(4) . The next objective is to develop an interactive redesign system that integrates 
optimization methods in order to have a flexible and efficient computational tool that 
can be used easily by design engineers. The interactive module is intended to create 
the missing link between BEM and CAD technologies and therefore it should constitute 
one of the key elements to computerize the shape optimization cycle. It has been found 
that the intrinsic definition of shape seems to suffer from a handicap that one doesn’t 
get the feel of the geometry of a curve. This problem can be eliminated by creating 
a multi-window facility with provisions for cartesian space, intrinsic space and various 
option buttons and information. The shape optimization cycle is composed of non- 
linear equation solver, numerical integration solver for mapping from intrinsic domain 
to cartesian domain, BEM module and a non-linear programming optimizer. 

(5) . Illustrative examples are to be considered and solved in order to show the 
capabilities of the proposed shape optimization methodology. However, emphasis will 
be focused on structural and mechanical applications. 

The organization of the dissertation can be outlined as follws: 

Chapter 1 introduces the relevent background information pertaining to the present 
research work. This includes a state-of-the-art review of optimization methods, an 
introduction of shape optimization, review of the relevant literature of intrinsic geometry 
and shape optimization problems, the objectives and scope of work, and the organization 
of the presentation of the dissertation. 

Chapter 2 deals with the shape synthesis using intrinsic geometry. The basic concepts 
of differential geometry are presented in this chapter. The intrinsic geometrical aspects 
of the LINear Curvature Element model of a curve are introduced. Some fundamental 
definitions such as curvature, torsion, arc length are discussed. The Serret-Frenet equa- 
tions which are the basis for curve synthesis are examined. For the case of planar curves, 
the differential equations relating to cartesian co-ordinates as a function of arc length 
are derived for the Serret-Frenet equations. The solution approach for synthesizing 2D 
curves are introduced. Intrinsic definition of the 3D surface is presented. Mapping pro- 
cedure from intrinsic to cartesian domain is discussed. Some illustrative examples are 
discussed to show the capabilities of intrinsic geometry based curve synthesis. Finally, 
some observation on 2D curve synthesis have been discussed. 
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Chapter 3 addresses the concepts of BEM analysis using intrinsic geometry. Geometry, 
load, and restraints specifications have been discussed. BEM analysis formulation for 
2D-elastostatic problems characterized by plane-stress, plane-strain, and axi-symmetric 
cases have been discussed. Sensitivity analysis aspects with respect to the SDVs are 
introduced. Finally, numerical example has been considered to show the numerical pro- 
cedure of stress-sensitivity with respect to the shape design variables and the graphical 
plot of the same has been shown. 

General approach of shape optimization methodology using intrinsic geometry and 
BEM has been addressed in Chapter 4. The integrating procedure between geometry, 
analysis and an exaustive search NLP technique has been discussed. Similarly the 
procedure of linking between geometry, analysis, and first-order optimization method 
has also been discussed. Issues of sensitivities of objective function with respect to the 
shape design varibles are addressed. The concepts of multi- window graphics environment 
to eliminate the handicap of intrinsic geometry has been discussed. Illustrative examples 
of some graphics environment have been presented. Finally, comaprision and efficiency 
of the proposed shape optimization methodology is made with respect to the existing 
methods. 

Case studies via Zero-order method are illustrated in Chapter 5. The illustrative 
examples using 2-LINCEs model include: a thick cylinder subjected to constant internal 
pressure and an elastic ring under diametral compression are presented. The illustrative 
examples using 3-LINCEs model include: a torque arm subjected to axial and transverse 
loading, a fillet under an axial load, and a ladle hook subjected to a tensile load. All 
these examples show how the intrinsic geometry based shape optimization can be used. 
An appreciable amount of reduction in the frontal area has been achieved depending 
upon the initial geometry of the component. Stress plots corresponding to the above 
examples are also presented. 

Case studies via first-order method have been presented in Chapter 6. As an illus- 
trative example the optimal shape of an elastic ring under diametral compression using 
3-LINCEs model is obtained. The problem formulation and results have been discussed. 
The problem formulation include establishing the constrained optimization problem by 
imposing suitable geometric and stress constraints and transforming to an unconstrained 
optimization problem using interior penalty function method. The DFP (Davidon- 
Fletcher-Powell) method has been used as an unconstrained optimization solver. 

Chapter 7 deals with the conclusions of the present work, which include a technical 
summary and the description concerninig the scope for future work. The use of intrinsic 
geometry a.s a geometric modeling tool in shape optimization has been evaluated. 
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Chapter 1 


INTRODUCTION 


1.1 An Overview of Optimization Methods 


The concept of optimiztion is intrinsically tied to humanity’s desire to excel. Though 
we may not consciously recognize it and though the optimiztion process takes different 
forms in different fields of endeavour, this drive to do better than before consumes 
much of our energy. Optimization is the process of maximizing a desired quantity or 
minimizing an undesired one. The process of mechanical design, or synthesis, may often 
be viewed as an optimization task. In optimization, the goal is to minimize an objective, 
such as weight or peak stress, subject to constraints on strength and performance. 
Optimiztion problems may be classified in to the following categories ([3], [41], [54]): 


(i) Standard or Size optimization. 

(ii) Combinatorial optimization. 

(iii) Shape optimization. 

(iv) Topology optimization. 


Standard or Size optimization : Size optimization relates to determining val- 
ues of design parameters such as thickness and cross-sectional areas and inertias. The 
mathematical optimization of a design assumes that a single measure of merit, objective 
function, or merit function exists and can be expressed formally in terms of the design 
parameters. The scalar merit function can be denoted by W^(X). Thus, the optimization 
problem can be stated as follows. 


Minimize W(X) 


(l.l) 
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subject to 


where 


^.(X) < 0, 


and hj{X.) = 0, 


z = 1, 2, ....m 

j = 1, 2, ....p 


X = 


Xi 

X2 

Xn , 


It should be noted that X is an n-dimensional vector called the design vector, VF(X) 
is called the objective function and 5 ',(X) and are respectively, the inequivality 

and the equivality constraints. The number of variables n and the number of constraints 
m and/or p need not be related in any way. The problem stated in Eqn. (1.1) is called 
a constrained optimization problem. 

If an optimization problem is simply stated as below: 

find yi = [xi,X 2 , - ■■ ,Xn)^which minimizes W(X.) (1-2) 

then such problems are called unconstrained optimization problems. 

It is likely that there may exist one or more values of X at which the function W will 
have a local minimum. Furthermore, it can be stated that the optimization problem 
is essentially to find an extremum value of the objective function. It can either be in 
the form of a minimum value of a function IT’(X) or it can also be in the form of a 
maximum value of — IT(X). 

There are many approaches to find the solution of an optimization. Basically, there 
are search methods that use the function values, VF, alone; gradient methods that 
use the function values and its first partial derivatives, and finally, methods that 

also require the second partial derivatives, The most common solution methods 

include ([3], [41], [54]) : 


(1) Search Methods. 

(2) Gradient Methods. 

(3) Second Derivative Methods. 

(4) Linear Programming. 

(5) Non-linear Programming. 
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Search methods include : exhaustive search, bisection search, grid search, random 
search, simplex search and golden search. Popular gradient methods are Steepest de- 
scent, Conjugate gradients, Newton’s method. The most popular methods include the 
method of Steepest descent, Conjugate directions, Powell’s direct search, method of least 
squares, and the Simplex search method. The search methods are easier to understand 
and are computationally inexpensive, since they do not require computing derivatives. 
Obviously, there are problems where utilizing derivative information will speed up the 
convergence to a local minimum. 

Combinatorial Optimization : In the standard optimization problems the decision 
variables are all continuous variables, that is they can assume all possible values within 
their range of variation. In combinatorial optimization, the problems considered are 
called discrete optimization problems in which some or all the decision variables are 
restricted to assume values within specified discrete sets. In combinatorial problems the 
optimum combination out of a possible set of feasible combinations has to be determined. 

If some or all the decision variables of a linear programming problem can assume 
only integer values then such a problem can be classified as a combinatorial optimiza- 
tion problem. Such problems are also called as integer linear programming or simply 
integer programming problems. The versatility of the integer programming model in 
several applications stems from the fact that in many practical problems, activities and 
resources, like machines, ships, and operators are indivisible. Also, most optimization 
problems of a combinatorial nature can be formulated as integer programs. 

Shape Optimization ; Shape optimization relates to determining the outline of a 
component and its dimensions such as, for example, the shape and size of a fillet or 
a hole. Shape optimization is a major consideration in many structural and mechani- 
cal engineering design problems. The problem of shape optimization for a component 
can be described as deciding the geometry of a curve or a surface which will maximize 
or minimize an objective function as well as satisfy a set of constraints. The exterior 
and interior boundary shapes of structures should be controlled and determined in the 
shape optimization process. Therefore, the parameters describing the changing bound- 
ary shapes of structures will be taken as design variables {referred to as Shape Design 
Variables). 

Shape optimization is more complex than sizing optimization ([3], [41], [54]). Since the 
shapes of structures or mechanical components are continuously changing in the design 
process, a number of points need careful consideration: describing changing boundary 
shape, how to maintain adequate finite element/boundary element mesh, how to enhance 
the accuracy of the sensitivity analysis, how to impose proper constraints, and how to 
utilize existing optimization methods to solve the shape optimization problems. 
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Topology Optimization : Topology optimization addresses the basic question of 
where material should and should not be. In various practical structural design prob- 
lems, decision of a layout of a structure or stiffeners in a given design space satisfying 
some design requirements is the most difficult task for design engineers since they must 
consider many possible conditions for safety, performance and so on. It is also very diffi- 
cult to solve the topology (layout) problem by any conventional structural optimization 
techinique because mathematical modeling of topology is too complex. To solve such 
problems, topology/layout optimization techniques are used based on various methods 
([19]). The topology optimization is applicable to actual design problems in automobile 
design such as frame and stiffener layout problems of a car body. 


1.2 Shape Optimization - An Introduction 


Engineering design is a methodology involving several stages starting with the state- 
ment of a need and ending with the final presentation of a design. Figure 1.1 shows the 
stages involved in the engineering design process. It can be observed that the design 
stages of synthesis, analysis and optimization constitute the core of any engineering de- 
sign methodology. The synthesis stage involves type synthesis, dimensional synthesis, 
material synthesis, synthesis of component tolerances as well as assembly clearances, 
and synthesis for manufacturability. The overall objective of dimensional synthesis is 
to design the geometry of physical objects. The geometrical aspects of the dimensional 
synthesis involve shape as well as size determination. The geometrical parameters could 
be described as scalar variables or as continuous functions. The diameter of a shaft is an 
example of a sizing design variable, the profile of a cam and a fillet of a structural mem- 
ber’s cross-section are examples of shape design functions. The shape design process 
can be carried out by co-operation between man and computer and can be visualized 
as shown in Figure 1.2. The shape design process can start by using a physical model 
and then convert it to a mathematical model. If necessary, the mathematical model can 
be modified, edited or updated. Moreover, one can create a physical model from this 
revised version of the mathematical model. 

The success of any shape optimization method depends on how the shape synthesis 
module (Figure 1.3) has been defined. The most frequently published approach is that 
of representing the shape in terms of either an array of control points or a spline like 
representation. This representation does not evaluate shape explicitly and does not 
indicate clearly how shape should be modified for improving a design. Although it 
is difficult to define, shape is still a basic property for design and manufacture. For 
example, the design of a cam or a fillet is the design of its geometry. The crucial 
decision is about the shape. Lord and Wilson (1984) defined shape as “Those aspects of 
geometrical form which have to do with the external aspect that an object presents to 
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the world”. Mathematical representation of a shape is difficult. Shape is like an overall 
attribute or characteristic. Unlike microscopic properties of a geometry like co-ordinates 
and slopes, shape is a macroscopic attribute. One can also find aspects of geometrical 
form that deal with extrinsic and intrinsic geometrical properties. The co-ordinates 
that represent the geometry are extrinsic shape properties and are analogous to mass 
and temperature and other thermodynamic and rigid body dynamic properties. Shape, 
on the other hand is independent of scale, position and orientation and can be viewed 
as an intrinsic property of geometry much like entropy is an intensive thermodynamic 
property. Examples of shape intrinsic properties are the curvature and torsion of a three- 
dimensional space curve. Figure 1.4 shows the role of shape in the engineering cycle. 
Shape is synthesized using a mathematical model. This is done in the intrinsic geometry 
domain. This representation of shape is then mapped into a cartesian domain, in which 
different analyses are performed based upon the particular application. The results are 
evaluated and updated till the optimum shape is found. Shape optimization is a field of 
research which is in rapid developement and is of increasing importance in engineering. 
For example, notches or holes of inappropriate shape, size, or orientation are known 
contributors to premature failure of mechancial components. Shape optimization can 
be used to determine the positions and sizes of notches or holes to minimize stress 
concentrations. 

Research into the optimization of structures dates back to the 17th century classi- 
cal treatise of Galileo, who investigated the optimum shapes of beams. In the 18th 
century, mathematicians posed certain classical shape optimization problems. For ex- 
ample, the “isoparametric problem” required finding the shape of a closed curve of fixed 
length that encloses the greatest area (answer: a circle). From 1960 through 1980, con- 
siderable research was done using classical mathematical approaches to determine the 
optimal shapes of columns, plates, arches, and shafts. The objective in these efforts 
was to maximize the buckling load (or lowest natural frequency) or to minimize the 
displacement in the structure for a given amount of material. 

In the last two decades, advances in FEM/BEM analysis using high-speed computers 
have shifted the focus to using numerical methods for shape optimization. Numerical 
methods are needed to optimize engineering components that generally have complex 
loading and geometry. Moreover, material reduction optimization needs simultaneous 
consideration of design and manufacturing. Shape optimization has been successfully 
applied to a variety of engineering problems in the automotive, aircraft and other indus- 
tries. A few examples are fillets, turbine disks and blades, pre-forms in forging, channels 
or ducts, stress minimum forms for elastic solids, and air-foil shapes. 

The method of shape optimization, proposed in this research work, consists of selecting 
a shape model in an intrinsic space, defining a set of shape design variables (SDVs) and 
then mapping into cartesian space. A shape model is represented as a set of LINear 
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Curvature Elements (LINCEs). The shape design variables are the values of curvature 
and/or arc length at some of the LINCEs. 


1.3 Literature Review 


The study of geometry is one of the most fascinating fields in the world of mathe- 
matics. Among the different branches of mathematics, geometry is probably the most 
appealing to our intuition. Historical reviews indicate that geometry was one of the 
fundamental branches of science developed by the Babylonians and the Egyptians. Eu- 
clid, one of the Greek Philosophores, established the foundation of what is known as 
the Euclidian Geometry. In the Ecuclidian Geometry, one studies the actual shape of 
objects. More accurately, one studies those properties of objects that are unchanged by 
rotation, translation and reflections of objects. 

Different branches of geometry have evolved in the past few centuries, projective ge- 
ometry and differential geometry are two examples. Both these fields are extensions of 
Euclidian Geometry. In projective geometry, one studies the way objects are seen in 
contrast to Euclidian Geometry in which one studies metric properties of real objects. 
Projective geometry has been of more interest to engineers than it has to mathemati- 
cians. It has become a tool for design engineers for conveying their ideas effectively. 
Differential geometry on the other hand deals with the geometry of curves and surfaces 
studied by means of differential calculus. Concepts of intrinsic geometry have been de- 
scribed in well known texts of differential geometry ([49], [13], [34]). In the literature on 
differential geometry, the Cornu spiral, known in American literature as a clothoid, but 
also called a transition spiral, or railway curve, since it is often used in rail track layout 
to increase or decrease curvature linearly is probably the widely illustrated example of 
a curve for intrinsic geometry. Intrinsic definition of a Cornu spiral is a linear curve in 
the curvature-arc length, k-s, plane. 

The developments in the area of computational geometry during the 70’s and 80’s gave 
an impetus for developing computational algorithms for many geometrical problems. In 
particular, curve and surface design became an active topic of research. One can now 
see a series of excellent texts which describe, in detail, how curves and surfaces can be 
designed using a variety of spline models ([17], [31], [34], [57]]). The major drawback of 
using a spline approach for shape synthesis is that the information about curvature and 
torsion is not explicitly available. This makes it difficult to provide a feedback from the 
analysis-optimization sequence to the shape synthesis algorithm. 

The works of Adams, Pal, Nutbourne, Remash Chandran, and Tavakkoli ([1],[37], 
[36], [42], [51], and [52]) form an interesting section of the literature reviewed here. The 
concepts of curve synthesis using linear curvature elements, termed as LINCEs, had been 
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developed by Nutbourne. This concept has been implemented in the form a detailed 
procedure by Adams (1975) and Nutbourne and Pal (1977) ([l], [37]). There are two 
difficulties in reviewing the works of Nutbourne, Adams and Pal. The first one being 
that they do not solve the boundary condition constraints of Serret-Frenet equations 
explicitly. This makes the computational procedure lengthy. The second aspect is that 
they have not shown how the method can be used for analysis-optimization process and 
hence sensitivity analysis. Remash Chandran (1982) ([42]) had used a spline curve to 
model the intrinsic shape definition. He had shown how a planar cam profile can be 
designed using their approach. It should be noted, however, that the method has not 
been extended for any shape optimization work. Tavakkoli and Dhande ([52]) worked 
on shape synthesis and optimization methodology using intrinsic geometry had shown 
only optimal configuration of Variable Geometry Truss design problems. It has not been 
extended for generalized shape optimization methodology involving FEM/BEM stress 
analysis and first-order optimization methods and hence sensitivity analysis. Although, 
the SDVs are the values of curvatures and/or arc lengths at some of the LINCEs during 
shape synthesis, the freedom and restrictions involved in choosing the number of arc 
length SDVs (i.e. s,s) as dependent design variables while solving the Serret-Frenet 
non-linear simultaneous equations and the robustness of numerical scheme have not 
been discussed. It is important to note that, out of three dependent shape design 
variables, the number of arc length (i.e s,s) variables can not be more than one. For 
example, in 2-LINCEs model, the design varibles will be Kq, sq, «i, ■Sii X’ti and S 2 . 
Out of these six SDVs by keeping Sq = 0, kq and si are chosen as free shape design 
varibles, and hence /Ci, si, and /ca will become dependent design variables which usually 
obtained by solving three non-linear Serret-Frenet geometric constraint equations, which 
is a successful case. Whereas Kq and ki are chosen as free SDVs, Si, K 2 , and S 2 will 
become dependent design varibles. Since the arc length variables are highly non-linear 
as compared to curvature variables, this choice will provide a difficult situation of not 
being able to converge to a solution. Such a behaviour of the Serret-Frenet equations 
has been verified using various non-linear equation solving schemes. The solution to 
the Serret-Frenet equations involve computation of Fresenal integrals. Unfortunately, 
the expressions of the Fresenal integrals are not available in closed-form. Hence, it is 
necessary to select a proper set of variables and also select an efficient as well as a robust 
numerical scheme for the solution of the Serret-Frenet equations. 

The finite element method has been extensively used for the analysis of mechanical 
components and in structural optimization during the last two decades, including suc- 
cessful application to shape optimal design of shafts, automobile components, aircraft 
and elastic structural components. In contrast, the boundary element method has only 
recently been applied to shape optimization of structural and mechanical components. 

The FEM application for shape optimization problems of structural and mechanical 
components has been successfully demonstrated, but at the same time it has some 
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drawbacks. It is often necessary to redefine new finite- element meshes as the geometry 
of the structure changes. Inaccurate evaluations of stresses at the boundary can be 
responsible for the evaluation of very inaccurate design sensitivity analysis, thus leading 
to a large number of optimization iterations or even unrealistic designs. These difficulties 
with the FEM can be partially overcome by using the BEM to discretize the structure. 
Results of the boundary element analysis of elasticity problems are more accurate than 
the corresponding solutions of the FEM models, and they are expected to yield improved 
design sensitivity information. Consequently, a smaller number of iterations is needed 
to find the optimum shape. In the last two decades large number of papers have been 
published in the evolution of the BEM for shape optimization of engineering design 
problems ([4], [9], [32], [33], [45], and [56]). 

The BEM is less versatile for structural analaysis than the FEM. Its applicability to 
shape optimization of structures is presently limited to elasticity problems subjected to 
static constraints. 

The merits of BEM against FEM can be summarized as follows: 

1. Less data preparation time. 

2. High resolution of stress. 

3. Less computer time and storage. 

4. Less unwanted information (since internal points in BE solutions are optional). 

5. Easily applicable to incompressible materials. The displacement based plane-strain 
FE formulation fails when Poisson’s ratio equals 0.5 exactly (i.e. the material is incom- 
pressible). The boundary element formulation, however handles these materials without 
any difficulty. Therfore, in problems involving rubber like materials the BEM is much 
more suitable than FEM. 

Perhaps a more severe weakness is the difficulty in modeling thin-shell structures using 
BEM. While BEM is capable of handling thickness ratios of 10:1 - and in some cases 
upto 25:1 - it generally does not produce accurate results for thin plates and shells. 

Traditionally, the optimization problem is cast as a minimization or maximization of 
an objective function subject to a set of constraints. A finite number of scalar variables 
are defined as design variables. However, if one needs to consider the geometry of 
a component as a design variable, then it no longer is a set of discrete variables. The 
geometry needs to be treated as a continuum. The problems of optimization dealing with 
continuum-type shape design variables are classified ais shape optimization problems. 
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An extensive literature is available on numerical methods for optimization of struc- 
tural and mechanical components whose shapes are defined by cross-section and thick- 
ness variables ([3], [41], [54]). Only limited literature has appeared in the area of shape 
optimization. The texts dedicated to this subject include: Pironneau ([39]), Haug, Choi 
and Komkov ([23]), Bennet and Botkin ([5]). The FEM has been applied extensively to 
shape optimal design of structures since 1973, while it is only in the last decade that 
BEM has been used in the field of shape optimization. 

Mota Soares, Rodrigues, Oliveira Faria and Haug ([32]) developed models for the 
shape optimal design of solid and hollow shafts, based on constant, linear and quadratic 
boundary elements and non-linear programming techniques. The design objective is to 
choose a shaft with a given area, which has maximal torsional stiffness. The models 
are much more efficient and robust than the corresponding finite element discretization 
since the sensitivity information is more accurate. 

Shape optimal design of shafts, based on the BEM, has also been developed by Bur- 
czynski and Adamczyk ([32]). Optimality conditions are generated and the Newton- 
Raphson method is used to solve a set of non-linear algebric equations. The examples 
show that the number of analyses required is smaller than a corresponding finite element 
discretization. 

Models for the shape optimal design of bidimensional elasticity based on the boundary 
element method and linear programming technique have been developed by Zychowski 
and Mizukami ([32]). The design objective is to minimize the area subject to displace- 
ment and geometrical constraints. 

Mota Soares, Rodrigues and Choi ([32]) have developed models for the shape optimal 
design of bidimensional structures based on linear and quadratic boundary elements 
and the non-linear programming technique. The design objective is to minimize the 
compliance having a constraint of constant area. Applications show that the BEM 
is more accurate and efficient than the corresponding FEM. For general shapes, the 
technique used to calculate the stresses with the boundary was not accurate. This 
problem has been overcome by Leal who also has developed an automatic technique for 
mesh refinement ([32]). 

Shape optimal design models for two and three-dimensional elasticity problems, based 
on BEM, have been developed by Burczynski and Adamczyk ([32]). The design objective 
is to maximize stiffness subjected to constant volume. The optimality conditions are 
derived for an optimal boundary. An iterative process, based on finite differences and 
Newton- Raphson method, is used to solve a set of non-linear algebraic equations, making 
it possible to determine the unknown optimal shape. 
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Eizadian and Trompette([32]) have also developed a model for shape optimal design 
of two-dimensional structures based on the BEM and NLP techniques. The design 
objective is to minimize the tangential stress subjected to geometrical constraints. The 
geometry is defined by linear and circular elements. Substructures are used to represent 
the fixed and moving boundaries. The multiplier method is employed to solve NLP 
problem. Several applications are presented including the shape design of connecting- 
rod and a rotor. Numerical instabilities are reported. 

The BEM has been applied to the shape optimization of air-foils and wings by Piron- 
neau ([39]). Furthermore, Fernand Ellyin ([18]) developed a methodology for shape 
optimization of intersecting pressure vessels. The methodology involves a profile of vari- 
able thickness which connects a spherical shell to a cylindrical one. The geometry of 
the mid-surface of the connecting shell of revolution is not known a priori, neither is 
the thickness variation. The profile and thickness is obtained based on minimum ma- 
terial volume and strength criterion. The optimum shape is obtained through a direct 
variation procedure of Rayleigh- Ritz method. 

Esping and Holm ([15]) published a research work on structural optimization using 
OASIS. OASIS is a code for structural optimization. This paper demonstrates the 
capabilities in shape optimization. The CAD system ALADDIN is used for shape de- 
scription and generation of finite element meshes. ALLADIN is a CAD system used 
to give meshes, loads, boundary conditions etc. to the finite element system, and to 
OASIS. It can also be used as a post-processor for stress, element errors etc. 

Hajela and Jih ([22]) have published their research work on boundary element method 
in Optimal Shape Design - An integrated approach. This paper describes an implemen- 
tation of BEM in optimal shape design. The computational advantages of the BEM over 
the more traditionally used FEM of analysis are examined in the context of the shape 
synthesis problem. An integrated analysis and optimization approach is also proposed 
for this problem. Numerical results are presented for shape design of torsional shafts 
and of structural components that can be characterized by a state of plane-stress, 

Sandgren ([44]) has published a paper on shape optimiztion - Creating a Useful Design 
Tool. In this paper several new areas of research are presented which seek to eliminate 
some of the barriers to achieving a useful design optimization tool. These areas include 
consideration of geometrical and topological optimization as well as cross-sectional op- 
timization, fundamental new approaches to combine analysis and optimization. Results 
for shape optimization using the boundary integral method with substructuring, zero- 
one integer programming for decision support and means of performing the analysis 
within the shape optimization formulations are presented. 
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Eschenaur ([14]) has published a paper on shape optimization of satellite tanks. The 
paper deals with the method of finding the optimal design of a thin walled satellite tank 
subject to constant internal pressure. 

Hou ([24]) developed techniques and applications of shape optimal design. This re- 
search objective of Hou is to present a unified formulation and efficient numerical tech- 
nique to perform shape optimum design. The numerical scheme consists of boundary 
parametrization, FEM and mathematical programming. 

Shape optimization of two-dimensional elastic structures with optimal thickness for 
fixed parts, has also been developed by Ding ([12]). Shape optimization of two-dimensional 
elastic sturctures is treated in this paper. B-spline function is used as a shape function 
and design element technique is adapted simultaneously. 

Boundary Integral Equation method for shape optimization of elastic structures has 
been published by Choi and Kwak ([9]). This work deals with a general method for shape 
design sensitivity analysis as applied to plane elasticity problems using a direct boundary 
integral formulation, material derivative concept and adjoint variable method. Among 
the several numerical implementations tested, the second order boundary elements with 
a cubic spline representation of the moving boundary have shown the best accurcy. 
Optimization as a philosophy of engineering structural design was reviewed by Schmit 
([51]). There is now a collection of excellent titles reviewing optimization methods in 
engineering design ([3], [41], and [54]). 

Shape optimization has attracted the attention of many researchers. Haftka and 
Grandhi have reviewed the work of several researchers in a review article ([21]). Van- 
derplaats has discussed the numerical problems associated with shape optimization. 
West and Sandgren ([56]) have used the variational approach to solve shape optimiza- 
tion problems. 

The concept of shape optimization has been applied to the design of many mechanical 
components. Some applications are mentioned here. 

As an application of optimization in automotive industry, Botkin([51]) used finite 
element analysis to design torque arm for an optimum shape which would minimize the 
weight subject to strength limitations. The optimum shape of rotating disks was studied 
by Bhavikatti and Ramakrishnan ([51]) using Non-Linear Programming (NLP) method. 
They used a 5th degree polynomial to define the shape of the cross-section. They also 
performed a stress analysis of the disk by FEM using isoparametric elements. They had 
satisfactory results which did not need any fairing process after the shape optimization. 

The optimum shpe of a gas turbine disk was obtained by Luchi, etal. ([30]). They used 
interactive procedure based on a method in which the disk profile was defined by spline 
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interpolation. Queue and Trompette ([40]) also worked on the design of an optimal 
shape of turbine blades. The weight of the rotating turbine blades was minimized in 
finding the optimal shapes of the cross-sections. They used polynomials to define the 
shape. 

West and Sandgren ([56]) investigated the optimum shape for a ladle hook using BEM. 
Seventh order B-spline curve with nine control points were used to define the shape. As a 
second approach, the shape optimization problem was formulated as a variational state- 
ment. A fillet problem was studied by Wu ([51]). The design objective is to find the 
optimum shape of a fillet in order to minimize the total area of the fillet. A fourth-order 
B-spline was used to represent the design patch. Wu obtained optimum shapes by spec- 
ifying different values of the permissible maximum stress. Mota Soares and Choi ([32]) 
have used boundary elements in shape optimal design of structures. They developed a 
technique for an automatic mesh refinement of the boundary elements. Tavakkoli and 
Dhande ([52]) have worked on the application of intrinsic geometry for shape synthesis 
and optimization to find a configurational design of a planar Variable Geometry Trusss 
(VGT). However, this work is not enhanced to integrate with FEM/BEM analysis and 
first-order optimization methods and hence with the sensitivity analysis. This work 
waa also not applied to the general application of engineering design problems involving 
weight or peak stress minimization. 

BESHOP [Boundary Element Shape Optimizer] - A Program Package for shape opti- 
mization in 2D-elasticity problems is developed by Farshi and Moradi ([16]). This work 
introduces a package which handles shape optimization in 2D-elasticity problems where 
BEM is used for analysis phase affording a large reduction in the problem dimensional- 
ity together with a variant of the method of center points for the optimization phase. 
Design variables in optimum shape determination consists of basic nodal co-ordinates 
for B-spline functions as the means to define smooth curves. Maximum stress and/or 
displacements at critical boundary points are considered as constraints in each cycle of 
the procedure. 

Belegundu ([6]) has worked on optimizing the shapes of mechanical elements. Accord- 
ing to his work, a set of basis shapes are generated using the pre- and post-processors. 
Then, an objective function such as weight or peak stress is minimized using a suitable 
optimization procedure. 

Design optimization using BEM has been published by Tafreshi and Fenner ([50]). 
This work addresses the integration between numerical optimization techniques and the 
BEM in order to optimize the shape of machine elements and structural components 
subjected to static loading. 
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Design Sensitivity Analysis using the BEM has been published by Tafreshi and Fenner 
([50]). This work deals with the application of boundary integral equations to design 
sensitivity analysis. Derivatives of displacement and stresses are calculated by implicit 
differentiation of the corresponding boundary element elasticity kernels. 

A review of the current state-of-the-art shows that the geometrical description of a 
component is taken into account either as a set of auxiliary points or as a set of points 
lying on the geometry of the curve or surface. Auxiliary points may not lie on the curve 
or surface. For example, the control polygon points of a Bezier or B-spline curve or 
surface are the auxiliary points of a particular curve or surface being designed. They, 
however, do define the shape of the curve or the surface. The review of the current 
literature shows that introducing the concept of intrinsic shape definition will definitely 
provide a better foundation for solving shape optimization problems. 

Among all the numerical techniques available for engineering analysis. Boundary Ele- 
ment Method is more suitable for shape optimization since one needs solutions only on 
the boundary. The shape optimization problems described above can also be formulated 
and solved using the shape optimization methodology based on intrinsic geometry pro- 
posed in this research work. The method of shape optimization, proposed here, consists 
of selecting a shape model, defining a set of intrinsic shape design variables (SDVs) and 
then evaluating cartesian co-ordinates of the curve. A shape model is conceived as a set 
of continuous piecewise linear segment of curvature (LINCEs), each segment defined as 
function of the arc length. The SDVs are the values of curvature and/or arc lengths 
at some of the Linear Curvature Elements. Then, mapping intrinsic space to cartesian 
space will be carried out. Thus, shape synthesis will be accomplished. The shape syn- 
thesis module thus obtained will be integrated through BEM analysis and an exhaustive 
search NLP (Non-Linear Programming) in one approach, and via first-order method of 
NLP in the other approach, to show how intrinsic geometry based shape definition can 
be used for engineering analysis and shape optimization. The proposed research also 
involves the method of calculating sensitivities of objectives like frontal-areas and the 
peak stress with respect to the intrinsic shape design variables. 

The optimal shapes of a thick cylinder subjected to constant internal pressure, an 
elastic ring under diametral compression, a torque arm subjected to axial and trans- 
verse loading, a fillet under an axial load, and a ladle hook subjected to a tensile load- 
ing are obtained as an application of this proposed shape optimization methodology. 
The complete proposed methodology has successfully been implemented as a CAD tool 
by developing multi-window concept for cartesian and intrinsic space using SUN-view 
graphics on SUN-workstation in a UNIX-environment and also on HP 9000/800-series 
workstation using STARBASE graphics on a UNIX-platform. In addition, it can be 
said that the proposed approach requires fewer shape design variables as compared to 
the methods where shape is represented using spline-like functions. 
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1.4 Objectives and Scope of Work 


The overall objective of this work is to develop a CAD (Computer Aided Design) 
approach to shape optimization using intrinsic geometry and Boundary Element Method 
and to show how this methodology can be applied to various structural and mechanical 
design problems. The specific objectives of this research work can be stated as follows. 

(1) . For the proposed methodology of shape optimization, it is required to develop 
and test a robust numerical scheme in order to solve Serret-Frenet equations. The Serret- 
Frenet equations play an important role in the curve synthesis. Since these are non-linear 
equations in intrinsic varibles, to be solved simultaneously in order to satisfy geometric 
constraints during the first stage of shape synthesis, proper caution must be taken in 
deciding free and dependent design variables. Suitable attention also must be given to 
select a stable numerical scheme and its implementation. The robust numerical scheme 
eliminates the designer intervention during the design process. This involves selection 
and development of a numerical integration module from amongst various numerical 
integration schemes such as Trapezoidal rule, Simpson’s rule and Gaussian quadrature 
etc. in order to develop a suitable scheme of shape synthesis and optimization. These 
numerical integration schemes will provide the necessary mapping from intrinsic domain 
to cartesian domain by solving the Serret-Frenet equations. 

(2) . It is also required to develop a generalized BEM analysis computer code based on 
existing sub-modules for 2D-elastostatic plane-stress, plane-strain, and axi-symmetric 
problems which take into account the geometrical definition of the boundary in the form 
of an intrinsic curve. Boundary Element Method for 2D-elastostatic problems will be 
applied as an analysis tool in the proposed shape optimization methodology and the 
problems considered are characterized by a state of either plane-stress or plane-strain. 

(3) . The key idea is to employ geometric modeling concepts of curve synthesis typical 
of the CAD technology in order to produce sensitivity analysis results. These sensitivi- 
ties with respect to intrinsic shape design variables can then be used by an optimizer to 
generate an improved design. Presently, analysis of engineering problems requires either 
FEM/BEM approach or either algebraic or numerical analysis methods. All these tech- 
niques are carried out in cartesian domain. The intrinsic geometry should provide the 
information about the cartesian geoemetry and any change in the shape made through 
the intrinsic design variables should update the cartesian geometry. In short, it is neces- 
sary to develop a method in which intrinsic geometry along with the cartesian geometry 
will be used for engineering analysis. To study and develop the numerical aspects of 
sensitivity of the curve or boundary of the component and hence the frontal-area of the 
cross-section and the magnitude of the peak stress (Maximum Von-Mises stress) with 
respect to the intrinsic shape design varibles. Determination of the sensitivity of the 
design to variations of the parameters is important in design and design optimization. 
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This seems particularly important since the intrinsic definition of the curve has been 
used to define the boundary of the component being designed in the present work. It is 
expected that changes in some of the intrinsic shape design variables will have a large 
impact on the resulting cartesian co-ordinates. 

(4) . The next objective is to develop an interactive redesign system that integrates 
optimization methods in order to have a flexible and eflhcient computational tool that 
can be used easily by design engineers. The interactive module is intended to create 
the missing link between BEM and CAD technologies and therefore it should constitute 
one of the key elements to computerize the shape optimization cycle. It has been found 
that the intrinsic definition of shape seems to suffer from a handicap that one doesn’t 
get the feel of the geometry of a curve. This problem can be eliminated by creating 
a multi-window facility with provisions for cartesian space, intrinsic space and various 
option buttons and information. The shape optimization cycle is composed of non- 
linear equation solver, numerical integration solver for mapping from intrinsic domain 
to cartesian domain, BEM module and a non-linear programming optimizer. 

(5) . Illustrative examples are to be considered and solved in order to show the 
capabilities of the proposed shape optimization methodology. However, emphasis will 
be focused on structural and mechanical applications. 


1.5 Organization 


The present research work is aimed at developing a CAD-methodology for shape 
optimization using intrinsic geometry and Boundary Element Method (BEM). The or- 
ganization of the dissertation can be outlined as follows: 

Chapter 1 introduces the relevent background information pertaining to the present 
research work. This includes a state-of-the-art review of optimization methods, an 
introduction of shape optimization, review of the relevant literature of intrinsic geometry 
and shape optimization problems, the objectives and scope of work, and the organization 
of the presentation of the dissertation. 

Chapter 2 deals with the shape synthesis using intrinsic geometry. The basic concepts 
of differential geometry are presented in this Chapter. The intrinsic geometrical aspects 
of the LINear Curvature Element model of a curve are introduced. Some fundamental 
definitions such as curvature, torsion, arc length are discussed. The Serret-Frenet equa- 
tions which are the basis for curve synthesis are examined. For the case of planar curves, 
the differential equations relating to cartesian co-ordinates as a function of arc length 
are derived for the Serret-Frenet equations. The solution approach for synthesizing 
2D-curves are introduced. Intrinsic definition of the 3D-surface is presented. Mapping 
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procedure from intrinsic to cartesian domain is discussed. Some illustrative examples are 
discussed to show the capabilities of intrinsic geometry based curve synthesis. Finally, 
some observation on 2D-curve synthesis have been discussed. 

Chapter 3 addresses the concepts of BEM analysis using intrinsic geometry. Geometry, 
load, and restraints specifications have been discussed. BEM analysis formulation for 
2D-elastostatic problems characterized by plane-stress, plane-strain, and axi-symmetric 
cases have been discussed. Sensitivity analysis aspects with respect to the SDVs are 
introduced. Finally, numerical example has been considered to show the numerical pro- 
cedure of stress-sensitivity with respect to the shape design variables and the graphical 
plot of the same has been shown. 

General approach of shape optimization methodology using intrinsic geometry and 
BEM has been addressed in Chapter 4. The integrating procedure between geometry, 
analysis and an exhaustive search NLP technique has been discussed. Similarly the 
procedure of linking between geometry, analysis, and first-order optimization method 
has also been discussed. Issues of sensitivities of objective function with respect to the 
shape design varibles are addressed. The concepts of multi-window graphics environment 
to eliminate the handicap of intrinsic geometry has been discussed. Illustrative examples 
of some graphics environment have been presented. Finally, comaprison and efficiency 
of the proposed shape optimization methodology is made with respect to the existing 
methods. 

Case studies via zero-order method are illustrated in Chapter 5. The illustrative 
examples using 2-LINCEs model include: a thick cylinder subjected to constant internal 
pressure and an elastic ring under diametral compression are presented. The illustrative 
examples using 3-LINCEs model include; a torque arm subjected to axial and transverse 
loading, a fillet under an axial load, and a ladle hook subjected to a tensile load are 
presented. All these examples show how the intrinsic geometry based shape optimization 
can be used. An appreciable amount of reduction in the frontal area has been achieved 
depending upon the initial geometry of the component. Stress plots corresponding to 
the above examples are also presented. 

Case studies via first-order method have been presented in Chapter 6. As an illus- 
trative example the optimal shape of an elastic ring under diametral compression using 
3-LINCEs model is obtained. The problem formulation and results have been discussed. 
The problem formulation include establishing the constrained optimization problem by 
imposing suitable geometric and stress constraints and transforming to an unconstrained 
optimization problem using interior penalty function method. The DFP (Davidon- 
Fletcher-Powell) method has been used as an unconstrained optimization solver. 
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Chapter 7 deals with the conclusions of the present work, which include a technical 
summary and the description concerning the scope for future work. The comparison 
and efficiency of the proposed shape optimization methodology has been made with the 
existing shape optimization procedures. Finally, the effectiveness of intrinsic geometry 
as a computer aided geometric design tool in shape optimization has been established. 



Chapter 2 


SHAPE SYNTHESIS USING 
INTRINSIC GEOMETRY 


2.1 Basic Concepts of Differential Geometry 


A few decades ago, researchers broke away from the main stream of research in compu- 
tational geometry from using polynomials of one sort or another and turned their ideas 
to the task of synthesizing curves and surfaces via differential geometry. Differen- 
tial geometry in its simplest form is the study of the shape of curves and surfaces using 
differential and integral calculus. Its origins lie in the writings of Euler and Monge 200 
years ago, and the subject has an immense literature spanninig many languages ([13], 
[26], and [34]). It provides a “study in the small or local study”, i.e. the properties 
of small piece of curves and surfaces. It can be explored that the intrinsic description 
of smooth curves as free curves in three-dimensional space. Arc length is denoted by 
s. At any point on a space curve there is infinity of normal vectors, but only one tan- 
gent vector. The tangent vector can be defined as the derivative of position vector 
t(s) = r (s) and that it has unit length. 

The derivative of the tangent vector defines both the principal normal vector n 
and the curvature k(s). The vector cross product of the tangent vector and the normal 
vector defines the binormal vector b(s). The tangent vector, the normal vector, and 
the binormal vector are three unit vectors all mutually orthogonal and define a curve 
frame f(s) (Figure 2.1) that moves along the curve as the arc length increases. The 
torsion r(s) to be introduced now should, roughly speaking, measure the magnitude 
and sense of deviation of a curve from the osculating plane in the neighbourhood of the 
corresponding point of the curve, or in other words, the rate of change of osculating 
plane (Figure 2.2). 

The three unit vectors t(s), n(s), b(s) form the space curve frame f(s) (Figure 2.1). 


23 



Figure 2.1 The Space Curve t-n-b Frame f(s). 
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f(.) 


rt(3) 

n(5) 


(2.1) 


The names given to the planes embraced by the vectors are as follows. The plane 
containing t and n is called the osculating plane. The plane containing n and b is 
called the normal plane. The plane containing b and t is called the rectifying plane. 
These planes are shown in Figure 2.2 ([13], [26], and [34]). 

Geometric elements are represented by use of either parametric or non-parametric 
equations. For a plane curve, an explicit non-parametric equation takes the form y = 
/(ar) which is not a convenient form to represent closed or multiple- valued curves. On 
the other hand, an implicit non-parametric equation takes the form /(x, y) = 0 in which 
multiple-valued curves are possible. In either case both these forms are axis dependent. 
Moreover, the slope at a point may become infinitely large depending on the orientation 
of the curve with respect to the co-ordinate axes. 

In parametric form, on the other hand, the slope values are computed in terms of the 
rate of change of the co-ordinates with respect to the parameter. For a plane curve, 
a set of two equations x — x(u), y = y{u), ^ u < Umax, with the parameter u, 

represent its parametric form. 

The formulation and application of Parametric Cubic (PC-curves), Bezier and B- 
spline curves have been thoroughly studied by Mortenson ([31]) and Rogers and Adams 
([43]). Neither one of the geometrical forms have the ability to control explicitly the 
intrinsic properties such as, curvature (^(s)) and torsion (t(s)) of a curve, which makes 
the shape synthesis using the intrinsic form more attractive. 

A space curve requires two intrinsic equations n = k{s) and r = r(s) in order to 
determine its unique shape. For the case of a 2-D curve, r = 0. Torsion r = 0 is a natural 
equation since it is a property of all 2-D curves. It suffices to claim that for any given 
function k(s), it is possible to solve the Serret-Frenet second-order equations for x(s) and 
y{s). However, only ramp type curvature functions which correspond to spirals in the 
cartesian co-ordinates have been selected in the shape synthesis methodology presented 
in this research work. The selection of multi-ramp functions to represent the curvatures 
in the Serret-Frenet equations provides the designer with more free parameters and 
hence more flexibility of controlling the shape. The effects of second and higher order 
curvature equations is not studied fully in order to explore the results of shape synthesis 
by means of these higher order equations ([42]). 

Serret-Frenet Equations 

Consider a curve C as shown in Figure 2.3. Let r be the position vector of a generic 
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point P and let s, the arc length of P from a reference point A, be the parameter 
describing the curve as r(s). The unit tangent vector, the curvature, the torsion, the 
normal and the binormal of the curve C at the point P are given as follows: 


‘=1 

(2.2) 

dt 

— — KXi 

ds 

(2.3) 

- = —/ct + rb 

(2.4) 

db 

ds 

(2.5) 

b = t X n 

(2.6) 


where t is the unit tangent vector, n is the unit normal vector, b is the unit binormal 
vector, K is the curvature, and r is the torsion. Equations (2.3), (2.4), and (2.5) are 
known as the formulae of Serret-Frenet([17], and [49]). These equations describe the 
location as well cis the orientation of a moving trihedron consisting of unit vectors t, n, 
b along the curve. These Serret-Frenet equations take a central part in the theory of 
curves. Furthermore, these equations can also be used as the basis for shape synthesis. 

If a curve is planar then r = 0 and the above equations can be written down as 
follows. 


a:(s) 

[y(^). 

(2.7) 

dv 

t = — and n • t = 0 
ds 

(2.8) 

dt 

— = /cn 
ds 

(2.9) 

dn 

— rr — /ct 
ds 

(2.10) 


In general, the torsion and curvature of a space curve can be directly related to the 
position vector r by the following equations ([17]). 
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where 


K = 


^ V 

ds ^ da^ 


TK 


2 _ I ^ ^ 


ds V ds"^ 


da^ / 



rfj dy dz 1^ 
ds ds ds J 


— - I £e ^ f 

^^2 ^ da"^ ds^ •* 


and 


d^r 

5^3 


d^x d^y 
ds^ ds^ 


is 

ds^ 


(2.11) 

( 2 , 12 ) 

(2.13) 

(2.14) 

(2.15) 


Solution Approach to 2-D Shape Synthesis 


The problem of finding the co-ordinates x and ?/ as a function of the arc length 
parameter s can be solved using Serret-Frenet equations (2.3) through (2.5). It should be 
remembered that the variation of the intrinsic property of curvature, k, gives flexibility 
to shape synthesis. 

To begin with, rewrite Eqn.(2.8) as follows: 

‘ = 1 = 1 '^ ( 2 - 16 ) 


Since n.t = 0, then 


n = 


ds 


(2.17) 


On the other hand : 


dt 

ds 


P d^x(s) 



ds^ 


Substituting Equations (2.17) and (2.18) in Eqn.(2.3) yields : 


(2.18) 
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Therefore, 




~M‘) 


d^y(s) 

= k(s) 

ds 

dxjs) 

(2.19) 

ds"^ 


^ ds ^ 


(Px 

ds^ 


= 0 

(2.20) 

(fy 

ds^ 

dx 

= 0 

(2.21) 


where the boundary conditions at the starting point A are assumed to be known i.e., 
at s = So, X = Xq, y = = ipQ. 


Assume that k{s) is defined as a piecewise linear function of s. Let 

/ t / 

T] = X + iy 


( 2 . 22 ) 


where rj' = ^, x' = g, and y' = 


The governinig Equations (2.20) and (2.21) can now be written as follows: 


(r/') - iK{s){ri') = 0 

(2.23) 

or 7^ = e J ^ ' = e ^ ^ 

(2.24) 

where / K(s')ds = 'tjj{s). 


Separating the real and imaginary parts, we get 


~ = cos[ip{s)] 

(2.25) 

^ = sm[^(3)] 

(2.26) 

Integrating these equations with respect to the arc length 
co-ordinates of the 2-D curve in terms of the parameter s. 

a yields the parametric 
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where 


x(s) = [ cos[il}{cr)]dcr + Xq 
y{s)= f sin[ip{a)]da + yo 

j So 



K{s)ds + '00 


(2.27) 

(2.28) 


(2.29) 


Equations (2.27) through (2.29) could be solved using numerical integration and nu- 
merical nonlinear equation solving for the cartesian co-ordinates x(s) and y{s) and the 
change in the tangent angle 0(s) along the 2-D curve once the curvature function k{s) 
is known. 


2.2 Intrinsic LINear Curvature Element (LINCE) Model of a curve 


Consider the problem of defining a curve passing through two points Po(a^O) Vo) and 
Pn{xniyn) in a two-dimensional space (Figure 2.4). It has also been assumed that the 
directions of tangents at Pq and Pn have been specified as 0o and 0n and the arc lengths 
from a reference point A as sq and Sn respectively. If « is the curvature at any point P, 
then it is assumed that the variation of k as a function of the arc length as the parameter 
K = k{s) has been specified. It can be seen that k{s) defines the shape of the curve. It 
should be remembered that the variation of the intrinsic property of curvature, /c, offers 
flexibility to shape synthesis. 

In this section, linear equations are used as shape models. Let the curvature k{s) be 
defined as a series of piecewise continuous linear functions with curvatures kq and at 
the initial and final points (Figure 2.5 and Figure 2 6). The intrinsic equation can be 
written as: 


(s) = 

//Cl - Ko^ 

j(s - So) ■+• Ko 

\ Si — So / 

(s) = 

1 i 

)(s - Si) + Kt 


So < s < Si 

■si < s < S2 (2.30) 


k(s) 



Sn— l) T ^n—l 
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Each linear segment corresponds to a spiral in cartesian co-ordinates. The area under 
the curvature equation represents the change in the tangent angles of the initial and 
final points of the 2-D curve. Therefore: 

^»-V.o= (2.31) 

Upon substitution of the curvature Equation (2.30) in Equation (2.29) and subse- 
quently in Equations (2.27) and (2.28) and considering the boundary conditions of the 
2-D curve results in the following equations. 


and 


Xn — 


i 


51 

sq 

52 


cos 


COS 


( — •®o<a') + Koa -f (7i 

I Si - So 2 

(— ^)(^ - -Sio-) + Kicr -h C2 

I S 2 ^ Si Z 


da 

da 


+ 

-b 


rsn 

Jsn^l 


COS 


(- )(— “ Sn-l<r) + Kn^ia + Cn 

L “ ^n-l ^ 


da 


(2.32) 


Vn - Vo 


fSl 

— / 3m ( 

JsQ L 


'30 

+ / sin 

Jsi 

+ 


Si - So 2 


{— — - Sia) -b Kia -b C 2 
. S2 — S\ 1 


da 

da 


-b / sm { )(-— - Sn-lO-) + -b C„ 

J»n-l L "Sn 'Sn-1 " 


da 


(2.33) 


where Ci through Cn are constants of integration. Note that Ci = ^o- The other 
constants can be exressed in terms of k,, s,, i = 0, ..., n and ‘ipo using the followig relation. 

^3-1 ~~ ■* J -2 ' ^ ' 

+ - /c,_is,-i, i = 2,3,..-,n 

Sj — Sj-i \ I J 

(2.33.a) 
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Figure 2.6 Shape Models : 

(a) 1-LINCE Model (R-Model n=l, 1=0) 

(b) 2-LINCE Model (R-R Model n=2, 1=0) 

(c) 3-LINCE Model (R-R-R Model n=3, 1=0) 

(d) 3-LINCE Model (R-R-R Model n=3, 1=1). 
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Examination of the right-hand side of Equations (2.31) through (2.33) indicates the 
existence of (2n-l-2) variables sq through Sn and kq through Since there are only three 
non-linear simultaneous equations to be satisfied, the left over free design parameters 
will be (2n — 1). Furthermore, if Sq is assumed to be zero, then the number of free 
design parameters reduces to {2n — 2). These (2n — 2) parameters can be termed as 
Shape Design Variables (SDVs). If a feasible set of these (2n — 2) SDVs are selected, it 
is possible to evaluate the shape as well as the objective function. In order to find the 
optimal shape, it will be necessary to search the feasible domain of these (2n — 2) shape 
design variables. 

The shape synthesis process can now be thought of as a two step process. The first 
step involves deciding the number n. Here, the selection of a shape model consisting of 
n linear curvature elements (LINCEs) in the intrinsic plane. Examples of shape models 
for n = 1,2 and 3 are shown in Figure 2.6. These models will be termed as R model, 
R-R model and R-R-R model for n = 1, 2, and 3 respectively. The second step of shape 
definition involves establishing the set of shape design variables for a given shape model. 

The total number of SDVs will be (2n — 2 — 1) where I could be the additional 
constraints relating the intrinsic variables ft(s and s[s. Notice that for n = 1 and I = 0 
the number of SDVs is zero and for n = 2, / = 0 and n = 3, / = 0, the number of SDVs 
would be two and four respectively. Furthermore consider a R-R-R model with n = 3. 
If an additional constraint of ki = «2 is imposed, then the total number of SDVs would 
be 3 since I will be 1. 

Once the total number of SDVs has been established, it is a matter of choice for the 
designer to select either all /c,s or s,s or a combination of some /c,s and some s,s as the 
shape design variables. It should be noted however, that there are two intricacies in 
choosing shape design variables. The first one is that, out of three dependent design 
variables which are usually obtained by solving three non-linear simultaneous equations, 
the number of .s,s should not exceed one design variable, since the Eqns. (2.32) and (2.33) 
are too non-linear in terms of s,s than k[s. Second difficulty is that, as the number n 
of LINCEs increases, the feasible domain of SDVs may become smaller and difficult to 
evaluate. Computationally stable schemes are being selected for solving Eqns. (2.31) 
through (2.33) simultaneously as well as for evaluating integrals of Eqns. (2.32) and 
(2.33). 


2.3 Intrinsic Definition of a Surface 

Consider a biparametric surface S shown in Figure 2.7. A generic point P on this 
surface can be described by the radius vector r(u, v) where the parameters u and v vary 
between 0 and 1. 
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r(u, v) 


x(u, v) 
y(u,v) 
z(u,v) 


0 < u,v < 1 


(2.34) 


The intrinsic properties of a surface patch can be described in terms of the unit normal 
vector at the point P, the tangent plane at the point P and the principal curvature at 
point P. The intrinsic quantities of a surface in the neighbourhood of a given point can 
also be described by the coefficients of the first and second fundamental forms denoted 
by Form I and Form II respectively ([31]). A brief description of these intrinsic properties 
of a surface are given below. 

The unit normal vector to the surface at point P is defined as 


n = 


r^ X r,, 
r„ X r„ 


(2.35) 


The tangent plane passing through the point P can be described by the following 
equation. 


(q - r)-(ru x ry) = 0 (2.36) 

where q is the radius vector of a generic point lying on the tangent plane. 

In order to define the curvature of a surface, it is necessary to define a specific curve 
lying on the surface and passing through the given point P. Consider a plane U perpen- 
dicular to the tangent plane and passing through the normal n. The curve of intersection 
between the surface S and the plane 11 is denoted as cr as shown in Figure 2.8. The 
curvature of the curve a at point P is defined as the normal curvature of the surface 
S at point P. 

If there is a curve C' obtained as the curve of intersection of an inclined plane 11, 
then the curvature of the curve a which is a curve of intersection between 11 and S is 
denoted by k. The curvature and k are related by the following equation. 


Kn = (/c ■ n)n (2.37) 

Now consider the effect of the variation of the curvature Kn 3S the orientation of the 
normal plane changes. Let a denote a reference angle in the tangent plane for describing 
the orientation of the plane 11. As a. varies from 0 to 27r, the curvature changes its 





Figure 2.8 Normal Curvature of a Biparametric Surface, 
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magnitude. The direction of the curvature vector Kn is along n for all values of a. It 
has Teen shown in the literature that within the range of 0 to 2r the curvature has a 
maximum and minimum value at values of a, say a* and a* + | ([49]). The extremum 
values of curvature are termed as the principal curvatures and the principal curvatures 
are denoted by k/ and /cjj. It can also be mentioned that the curvature of a surface is a 
tensor and analogous to a quantity such as the stress at a point. In short, the properties 
such a5 the unit normal vector n, the principal curvatures /c/, kji and the angle a* 
uniquely characterize the intrinsic geometry of the surface S at point P. 

The intrinsic properties can also be described by means of coefficients E, F,G, and L, 
M, N where : 


and 


E = Vu-Tu 


F = r„-rt, 

(2.38) 

G = ri,-r„ 


L = ruu-n 


M = ruu-n 

(2.39) 

N = r,;„-n 



These coefficients are used to express the well known first and second fundamental 
forms as follows: 


FormI = Edu^ + 2Fdudv + Gdv"^ 
For ml I = Ldu^ + 2Mdudv + Ndv"^ 


(2.40) 


Furthermore, it has been seen that the principal curvatures are the roots of the following 


equation. 


[EG - F^y - {EN + GL- 2 FM)k + {LN - M') = 0 (2.41) 

where E, F, G and L, M, N are the coefficients described in Eqns. (2.38) and (2.39). 


The foregoing description about the intrinsic properties of a surface needs to be taken 
in to account for designing or synthesizing the geometry using the intrinsic properties. 


2.4 Intrinsic-Cartesian Mapping 


Two numerical schemes have been used in the present work and these are described 
here as a ready reference. These robust numerical schemes are used for the numerical 
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integration and solving non-linear simultaneous equations (2.31), (2.32), and (2.33) in 
order to carry out intrinsic-cartesian mapping. The numerical schemes are as follows. 

1. A Numerical Integration Scheme based on Simpson’s | Rule. 

2. A Nonlinear Equation Solver based on Newton-Raphson Iterative Scheme. 

Numerical Integration Based On Simpson’s Rule 

Let y = f{x). In order to evaluate f(x)dx, it is required to follow the procedure 
described below. 

Divide the interval from a to b into 2n equal parts each of width h, and denote the 
ordinates 


f{a),f{a + h),f{a + 2h), ,/{a+ (2n - l)h},f{b) 


as 


Then 


2/0?yi)J/2) 1 2/2n-"l) J/2n • 


1 

J f{x)dx = -h[yo -f ?/ 2 n + 2 (t /2 + yA+ + y2n-2) 

+ 4(T/i+l/3-f + y2r!-l)]- 

Consider, a typical integral term as given in Equation (2.32). 

Let 


(2.42) 


/ = / cos[(-^^^^ — —){^-SQ(j) + Koa-\-CA]da 

J sq *3 1 5o " 

This integral can be evaluated using Equation (2.42). The pseudo-code of a procedure 
required for this purpose is given below. The Same procedure can be used to evaluate 
other integral terms of Equation (2.32) and (2.33). 


double int_cos(/co, Ki, So, Si, ni, c) 

double « 0 ) ■So, si; 
int n\] 
double +c; 
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{ 

double a,b, sum, temp, hi] 
int t] 


a = So] 


b = 5i; 


hi = 


(2*ni) ’ 


sum = 0; 


sum = sum + cos[[{ki — Ko)/{si — so)) * ((o * u/2) — so * o) 
+ «o*o + +(c + l)) 

+ COs(((ki — K.o)/{Sl — 5o)) * {{b * 6/2) — So * 6) 
4- Ko * ^ + 1)); 

for(i = 1; i <— (2 * ni — 1); *+ = 2) 

{ 

temp = a + i * hi] 


sum = sum + 4 * (cos(((ki ~ kq)/ ( si - Sq)) * {(temp * temp(2) sq * temp) 
+ Ko * temp + *{c + 1)))] 


for(i = 2; i <= {2*ni — 2); i+ — 2) 

{ 

temp = a + i * hi] 

sum = sum + 2* (cos(((«i - ko)/{si - so)) * ({temp * templ2) - sq * temp) 
+ Ko * temp + *(c 1)))] 


} 

sum = sum*(hi/3)] 
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return {sum)] 

} 


Nonlinear Equation Solver Based On Newton-Raphson Iterative Scheme 

It is frequently required to determine the value of x at which a function f{x) is 
zero. This value is called as the root of the function /. A commonly used method for 
determining the root is the Newton-Raphson scheme. It estimates the new value of the 
root by using the local tangent direction of the function evaluated at the assumed value 
of the root. Improvement in the estimates is carried out iteratively. This is shown in 
Figure 2.9. 



Let P be the point at which / is zero; let xj be the initial guess for the root, fi = /(xi) 
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and and ^ a.t x = xi. The improved value of the root ( 2 : 2 ) is obtained from the 
formula, 

tanO = f[ = -r -^ — r 

{xi-a:2) 

Hence, 

X2 = Xx-[^). (2.43) 

h 

Iterations continue till / calculated at is reasonably small. The procedure converges 
only if the initial guess for x at P satisfies the condition, /( 0. The above procedure 

converges rapidly as the root is approached. 

If /(x) is known as an algebraic expression, /'(x) can be determined in closed form. 
If this expression is not available or is too complex being differentiated, the derivative 
is calculated as. 


where 


8 ~ a small value 


The roots of a pair of simultaneous equations. 


x,y) = 0 
x,y) = 0 

is determined as follows. Moving along the tangent direction in the one variable problem 
is equivalent to moving along a tangent plane (or surface) in a many variable problem. 
On this plane the tangent is represented by the total differential of the function. Using 
the rules of partial differentials we get. 



where 


d/ « A/ = /lAx + fyAy 

dg^ Ag = g^Ax + gyAy 


(2.45) 


f = ^ f 

8x1 Jy 


etc. 


If X and y are the present estimates of the roots, we need to find the corrections Ax 
and Ay such that x + Ax, y + Ay make / and g as zero. This gives Af = /(x + Ax,i/ + 
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Ay) - f{x, y) ~ -/{i, y) and Ay = g{x + Ax, y + Ay) - y(x, y) = -y(x, y). The two 
equations given above can be solved for Ax and Ay to yield, 


Ax = 
Ay = 


(J9u~9fv) 

det 

{gfx-fgx) 

det 


(2.46) 


where 


dei — {9xfy fx9y)~ 


Corrections Ax and Ay must -be repeatedly applied to x and y till convergence is 
attained. This requires the magnitudes of x and y to become small or the values of / 
and g to be close to zero. 

The above procedure can be generalized to n-simultaneous equations /i = /2 = • ■ • = 
/n = 0. The corrections to the roots xi, X 2 , • • • , x„ are calculated by inverting the matrix 
equation. 

The A matrix is given as A = [/jj], i,j = l,2,-'-,n where /,j = Matrix 

inversion can be carried out using Gaussian Elimination method. As before, iterations 
continue till Ax,’ is small and/or /,s are close to zero. 

It is very important to note that, lot of care must be taken while selecting the de- 
sign variables during Intrinsic-Cartesian mapping. It is due to the fact that the Eqn. 
(2.32) and (2.33) are too nonlinear in terms of s,s compared to k,s. It has been tested 
through various non-linear equation solvers and Generalized Leibnitz rule ([35]) for par- 
tial derivatives of Equns. (2.31), (2.32), and (2.33) and concluded that, out of three 
dependent design variables which would be obtained by solving above three equations 
simultaneously, the arc length variables (s[s) can not be more than one. In short, for 
example in a 2-LINCEs model select ki, K 2 , and as dependent design variables and 
hence kq, Si would be free shape design variables and So = 0 is always assumed. Instead, 
if Ki, si, and S 2 are chosen as dependent design variables, then it is impossible to solve 
equations (2.32), (2.33), and (2.31) simultaneously and the numerical process will never 
converge at all. 

In order to solve three non-linear Equations (2.31), (2.32) and (2.33), the Newton- 
Raphson method applied to 2-LINCEs model intrinsic-cartesian mapping can be seen 
as a detailed procedure described below. 

Consider a 2-LINCEs model with the following shape design variables (SDVs), ko, Sq, 
Ki, si, K 2 , and S 2 . Out of these (2n -|- 2) SDVs, it has been assumed that sq = 0, and 
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kq and are selected as free SDVs, hence ki, K 2 , and S 2 would be dependent design 
variables. Then by applying the Newton-Raphson numerical scheme as explained above, 
we get the following. 

Let 


f = fii^i,S 2 , K 2 ) represent Eqn. {2.32) 

9 = g{i^iiS 2 i K 2 ) represent Eqn.{2.33) 
and w = w{ki,S 2 , K 2 ) represent Eqn. {2.31). 

Therefore, 

/«i A/ci + /S2AS2 + /K2A/C2 = 0 

+g,^As2+gK2AK2 = 0 

U;«, AKi + As2 + 10kjAk2 = 0 

where 


(2.47) 

(2.48) 

(2.49) 

(2.50) 

(2.51) 

(2.52) 


f =-AL f 


df 


9«2 


w 


«1 


Ski ’ 


3w 


dw 

dK2 


For a 2-LINCEs model Eqns. (2.47) through (2.49) can be written as 


f = X 2 ~Xq - / cos[( 

•/ SO 

- / cos[( 

''*1 


..K2 — Ki 


)(- - SqCt) + Kocr + CildcT 
Si — So 2 
2 

)(- - Sict) + Kja + C'2]dcr = 0 


S2-Si^^2 


(2.53) 


9 = 92-90 - / sin[{ )(--soo-) + A:o^r + Ci]c(cr 

J 3q tSi — 5o ^ 

2 

- r sm[( - - ~ - - )(- - si<t) + Kia + C 2 ]d(T = 0 (2.54) 

^2 ^ Si 2 

«. = V -2 - V-o - - »o) - -s:) = 0 (2.55) 

The Eqns. (2.50) through (2.52) can be represented in matrix form as follows : 


/ki /s 2 f K2 


' Aki 


■ -/(Ki,S2,«2) ' 

9f^l 9 s2 9k2 


As2 


-^(/Ci,S2,K2) 

_ Wsg, ^ 




_ -U;(«1,52,K2) _ 


(2.56) 
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A^i = (—Aif — A^g — Azw)ldetA 
As 2 = {—Bif — B2g — Bzw) I detA 
A/C2 = {—Cif — C2g — Czw)ldetA 


where 


Ai = gs^w^^ -gK2Ws2 
Bl = -gK^Wn2 

Cl gK\^S2 9 s2^K.I 

A 2 = -fs2WK2 + Sk2W,^ 
B2 = + fK,WK2 

C 2 = -f^iW,2 + Ss2W^2 
A3 fs29K2 fK29s2 

B 3 ~ fril9K2 

C3 — “■ /s25«1 


detA — fniAi + /*2 -Bi + /kjC*! 


The new values of design variable k\s and/or s|s would be updated as follows. 


^ + Aki 

52 '4= S 2 T As2 
K2 4= /C2 + Ak2' 

The iterative process is repeated till the convergence is obtained and/or function 
values are close to zero. 

Algorithm 

Since the proposed methodology is a CAD approach to shape optimization, the shape 
synthesis and hence optimization module must be accompanied by an interactive graph- 
ics environment. This is an important issue since the intrinsic profile doesn’t give the 
real feel of the geometry to the designer. In order to eliminate this handicap of the 
intrinsic geometry, the first step is a creation of multi-window interactive graphics en- 
vironment that will include one window for the cartesian space, one window for the 
intrinsic space, the third one to have menu section buttons, and the fourth one to show 
the necessary written information or feed back within the graphics environment. To 
create this nature of graphics environment, computer program written in C has been 
developed based on SUN-view graphics on SUN-workstation with UNIX-platform. The 
next step in shape synthesis is to select the feasible shape model such as, 2-LINCEs, 
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3-LINCEs model etc. based on the geometry of the problem. Once, the shape model has 
been decided, the design variables can be specified a.s the input. A computer program 
has been developed in C using Simpson’s | rule, Newton-Raphson non-linear equation 
solver as described above in order to solve Eqns. (2.31), (2.32), and (2.33) simultane- 
ously and ensuring free and dependent design variables are selected properly. Thus, once 
all the shape design variables are established for a particular LINCE model, the next 
step will be to map from intrinsic to cartesian space. Thus, the intrinsic geometry based 
shape synthesis would be accomplished. This shape synthesis module would be used as 
a subroutine in any CAD based shape optimization cycle. The input being specified to 
this shape synthesis module in addition to the intrinsic variables k\s and/or s(s would 
be, the initial and final co-ordinates (i.e. {xo,yo) and [xn,yn) ) of a curve and the initial 
and final tangent angles (i.e. ■00 and 0„) of the corresponding points. The final step of 
this intrinsic-cartesian mapping would be the displaying of the cartesian and intrinsic 
profiles of a particular component or structure. The effect of change in intrinsic design 
variables can be seen in cartesian space immediately through the multi-window graphics 
environment. The data files that are created for the display purpose have been devel- 
oped as script files, so that they can be compatible with any graphics or plotting aids 
such as AutoCAD etc. 

2.5 Illustrative Examples 


When shape synthesis is carried out using intrinsic geometry, the designer has to select 
(i) appropriate end conditions, (ii) a suitable shape model, and (iii) an appropriate set 
of free and dependent shape design varibles. This has been illustrated by means of 

1- LINCE, 2-LINCEs, and 3-LINCEs model examples. 

In Example 2.1, it is required to define a curve passing through points (0,0) and 
(12,10). The tangent angle at the first point (0o) is zero. The tangent angle at the 
end point (0„) is taken to be 120 degrees. In this case, 1-LINCE model (R-model) is 
selected as the shape model. The variation in intrinsic variables can be seen as dotted 
curves in both intrinsic and cartesian space. The result has been shown in Figure 2.10. 

The selection of different free parameters can be illustrated with an example 2.2. A 

2- LINCEs model (R-R model) is selected as the shape model in this example. The shape 

design variables for this model are kqi -Si S'Hd S 2 . The solid curve is designed 

by choosing Kq and si as free SDVs and ki,S 2 , and «2 dependent SDVs. The values 
chosen for Kq and si are 0.6 and 12 respectively. The co-ordinates of end points are 
given as (0,0) and (12,10) and the tangent angles 0o = 0; 02 = -70 degrees. The Figure 
2.11 shows the result of this curve design. 

Example 2.3 is shown in Figure 2.12. The curvature model is selected being a 3- 
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LINCEs model (R-R-R model). The initial co-ordiates are at (0,0) with zero degree 
tangent angle (^o = 0 ) and the final co-ordinates are ( 12 , 10 ) corresponding to the final 
tangent angle of zero degree. In this curve design problem, the SDVs include «o, «i, 
Si, K 2 , sj, K 3 and S 3 . Out of these SDVs, « 2 j ^3 and S 3 are chosen as dependent design 
variables and hence kq, ki and Si becomes free SDVs. 


2.6 Some Observations 


It should be noted that when a single ramp function ( 1 -LINCE model) is selected 
as a shape model, no free shape design variables exists. In such a case. Equations 
(2.31) through (2.33) result in one unique answer. It has been found out that out of 
three dependent design variables, which are usually obtained by solving Serret-Frenet 
Equations (2.31) through (2.33) simultaneously, the number of arc length varibles (i.e. 
s,s) can not be more than one. For example, in a 2 -LINCEs model the selection of /Ci, 
/C 2 and S 2 as dependent design variables would work perfectly. Instead, K 2 , si and S 2 or 
Ki, si, S 2 or Kq, Si and S 2 combination of dependent design variables would not be able 
to give converged solution by using any non-linear equation solver. It is owing to the 
fact that the Equations (2.32) and (2.33) are too non-linear in terms of s,s compared to 
«,s. The situation will be similar in case of 3-LINCEs or n-LINCEs model. 

The numerical schemes which have been selected for this shape synthesis are found to 
be more robust and stable. This eliminates the selection of specific input for shape syn- 
thesis as well as the intervention of the designer during the design process. Any mechan- 
ical and structural engineering shape optimization problem requires a shape synthesis 
module. Traditionally, this has been accomplished by using parametric spline functions. 
The co-ordinates of the control polygon vertices of the spline functions constitute the 
shape design variables. Such an approach invariably leads to a large number of design 
variables. In contrast, the intrinsic shape definition clearly indicates that the number of 
shape design variables is definitely less than the number of shape design variables in the 
parametric approach. This intrinisic geometry based shape synthesis module is an asset 
to the shape optimization cycle and will become an effective computer aided geometric 
modeling tool. 


central liBRAW 

t I. T., KAWPUB 

A. rnSommi^t^ 
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Since, there are 2N unknowns, it need 2N equations to solve the problem. Looking 
at it in a more general way, it requires N sets of equations where each set contains two 
equations. For example; a force (or traction) placed at node 1. Using the fundamental 
solution, it is possible to calculate displacements and tractions at every node from 1 
to N. This yields a first set of linear equations (i.e. rows 1 and 2 of the matrix). To 
produce a second set of linear equations (i.e. rows 3 and 4 of the matrix), place the 
force at node2 and repeat the use of the fundamental solution to calculate the variables 
at the other nodes. This operation is repeated until the force is placed at the last node 
N, which will give the final set of equations (i.e. rows 2N — 1 and 2N of the matrix). 
Therefore, we will end up with 2N equations and 2N unknowns which produces a unique 
solution. 

The data required to fully define a given problem must include the following: 


(a) Data such as the number of Gaussian points, two-dimensional (plane-stress or 
plane-strain or axi-symmetric) geometry and material properties. 

(b) Mesh data such as the co-ordinates (in cartesian space) of the nodal points and the 
nodes on each element. 

(c) Boundary conditions such as prescribed displacements and/or stresses. 

(d) Co-ordinates of any internal points. 

The description of the above mentioned data is shown below: 


ANALYSIS DATA 

Geometry; This defines the type of geometry : either 

1 (Axisymmetric) 

2 (2-D plane stress) or 

3 (2-D plane strain) 


Gaussian: This difines the number of Gaussian integration points used for numerical 
integration. Either four or six Gaussian points are allowed. The use of four Gaussian 
points is sufficient for good accuracy. However, for thin structures better accuracy is 
obtained with a larger number of Gaussian points. 

E : Young’s modulus 

i' : Poisson’s ratio 


NODAL CO-ORDINATES (Figure 3.2) 
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Node 


: Node number 


X : X - co-ordinate 

y : y - co-ordinate 


ELEMENT DETAILS 


Element 
Node 1 
Node 2 
Node 3 


: Element number 
: First node of the element 
: Second node of the element 
: Third node of the element 


DISPLACEMENT BOUNDARY CONDITIONS 


Element : Element number 

Direction : The direction of the prescribed displacement : either 

1 (x - direction) or 

2 (y - direction) 


Disp 1 : Displacement of the first node of this element 

Disp 2 ; Displacement of the second node of this element 

Disp 3 : Displacement of the third node of this element 


STRESS BOUNDARY CONDITIONS 


( cTxx, <^yyi o'xy ) 1 : Prescribed stress on the first node of this element 

( cTxxi CTyy, cTxy ) 2 : Prescribed stress on the second node of this element 

( cTxxi ^yyi <^xy ) 3 : Prescribed stress on the third node of this element 


(Note: Tensile stress is considered positive, while compressive stress such as pressure 
is considered negative. If a normal pressure, P, is prescribed on a curved surface, simply 
input CTxx = (Tyy = -P and o^y =0.) 

INTERNAL POINTS 

Point : Internal point number 

Xint : X - co-ordinate of the internal point 

yint '• y - co-ordinate of the internal point 
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3.2 BEM Analysis 


Betti’s reciprocal work theorem and the Somigliana’s identity ([4]) for the displace- 
ments can be used to derive an integral equation for the displacements at a interior point 
p due to tractions and displacements on the surface at a boundary point Q (as discussed 
in Section 3.1). In the absence of body forces, the Boundary Integral Equation (BIE) 
can be written as follows: 


Ux{p) 

Mp) 


. f \ Txx{p,Q) Txy{PiQ) 

Jr [ Tyx{p,Q) Tyy(p,Q) 

_ f Uxx{p,Q) UxyiPtQ) 
Jr _ Uyx{piQ^ ^yyijPiQ^ 


Uy(^Q) 


dr{Q) 


tx{Q) 

iyiQ) 


dr{Q) 


(3.1) 


Using tensor notation, the above two BIE equations can be expressed as follws; 

Mp) + 1 TUp, Q)u,{Q)dr{Q) = j^u„{p, Q)t,{Q)iW) (3-2) 

Before attempting the numerical implementation of BIE, it is necessary to examine 
its physical significance. The fundamental solution (represented by the kernels 17, and 
Tij) simply works out the “ influence ” of a concentrated point force at a given point p 
on any other point Q (except when p coincides with Q then it becomes singular). This 
fundamental solution is valid for any elastostatic problem regardless of its geometry. The 
BIE of Eqn. (3.1) effectively relates two distinct problems; the first is the fundamental 
solution of a point force (the kernels Utj and T,j) and the second is the problem of 
solving the unknown variables u, and tj. 

Numerical Implementation 


In order to solve the BIE numerically, the steps required are as follows: 


STEP 1 : Division of the Boundary into Elements 

The boundary of the solution domain must be divided into a number of connected 
elements. Over each element, the variation of the geometry and the variables (dis- 
placements and tractions) must be described. These variations can be constant, linear, 
quadratic, cubic or higher order. Furthermore, it is possible to allow the geometry vari- 
ation to be different from the variation of the variables (for example a linear geometry 
description but with a quadratic variation of displacement). The literature review on 
BEM and FEM shows that, in most stress analysis applications, isoparametric quadratic 
elements provide the best compromise between accuracy and efficiency. Isoparametric 
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elements are elements that use the same order of variation for both the geometry and 
unknown variables. The elements which are used in the stress analysis of this research 
are of isoparametric quadtaric elements, but the change in the order of variation can be 
made by prescribing a new ‘shape function’. 

Shape functions are the most convenient way of describing the behaviour of any 
element because they use number of points on each element ( called ‘nodes’ or ‘nodal 
points’ ) where the variable value is given. Therefore, for a quadratic element, it require 
three nodes on each element: one at the mid point and one at either end. It can be 
defined a new co-ordinate system that is local to the element using a variable, with 
its origin at the mid point of the element and values of -1 to +1 at the end node (Figure 
3.2). 



Therefore, the geometry of an element can be described by the co-ordinates of its 
three nodes and the shape function as follows: 
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3 

+ -^2(0*2 4" ^^{ 0^3 

C=1 

3 

yiO = E ^ciOVc - NriOvi + N2{0y2 + N^iOys (3.3) 

C=1 

where the shape functions NdO are quadratic functions that must satisfy two condi- 
tions; 

!• ^c{0 = 1 at the node c and 
2. Nc{^) = 0 at the other two nodes. 


The explicit expressions for the shape functionsare as follows: 

m) = 

m) = (i+f)(i-e) (3.4) 

m) = f(i+{) 

Similarly, because the elements are isoparametric, the shape functions are used for 
the solution variables, follows: 


MO = 

C=1 

= 

MO = ZMOiM^ 

= MOMi + MOiOyh + mXMs (3.5) 

MO = tmXMc 

= k(0(f.), + mx‘.)2 + iV3(0(i,)3 

MO = LMOi/.X 

= rnOMi + + mXM3 
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The explicit expressions for the shape functions in terms of ^ are given in Equns. (3.4). 
STEP 2 : Numerical Integration of the Kernels 

The numerical integration is now performed over each element using the local co- 
ordinate ( rather than the boundary curve parameter T. 

Because of the transformation of the variable from the boundary curve T to the local 
co-ordinate if, the Jacobian of transformation J must be calculated as follows: 




di 


di 


+ 


\dy[i) 


d( 


(3.6) 


In order to determine the components of the unit outward normal, the unit tangential 
vector, m can be defined as follows: 


Unit Vector m = 7—763; + 7~7fiy 
|m| |m| 


where the magnitude of the vector m is given by 


ml = + (mj/)2 = 



and ex, By are the unit vectors along the x and y axes respectively. 


(3.7) 


(3.8) 


The magnitude of the vector m is equal to the Jacobian J(f), defined in Eqn. (3.6), 
Therefore, the components of the unit tangential vector can be written as follows: 


rUx = 
TUy = 


df J 

1 r dy(f) i 

df J 


(3.9) 


Consider now the unit vector in the z-direction, (normal to the two-dimensional xy 
plane). The unit normal vector is, therefore, equal to the vector product (cross product) 
of the vectors m and as follows: 


n = m X 02 (3.10) 

where 655 is the unit vector along the z-axis. The component of the unit outward normal 
(of unit length) are given by 
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Ur 


n„ 


1 


dy{i) 


-1 rdx{() 


d^ 

The differentials of the co-ordinates a;(^) and y[^) with respect to are given by 


dx{0 dN^iO , dN,{0 , dNsiO 

-If = -IT""' ^ ~ir"" + 

Mi) <im) . dN,{0 , dN,{i) 

— yi + —if—yi + 


(3.11) 


di d( d(, di 

where the differentials of the shape functions are given by 


(3.12) 


dm) 1 

di ^ 2 

di 

mio 1 

2 


(3.13) 


The BIE of Eqn.(3.1) can now be written as in terms of the local co-ordinate ^ as 
follows; 


+ 


Cxxij?^ ^xyi^P^ 

■ Wx(p) ■ 

. Cyx{p) Cyyip) . 

. ^y{p) . 


\ 'd^XxiPiQ) Txy{p,Q) 

[ TUp.Q) T„(p,Q) 

vMp,Q) 

lr,.(p,Q) u,.{p,Q) 


N.imodi) 

m)3(i)di) 


■ Ux(Q) 

. Uy{Q) . 

■ ix{Q) ■ 
L ty{Q) J 


(3.14) 


where M is the total number of elements. The integrals in the above equations can be 
placed in new functions [A] and [B] as follows; 


Cxxip) Cxyip) 
^yx{p) ^yyi.P) J 


Ux{p) 

«v(p) 


-^xy 


■ Ux{Q) 

. Ayx Ayy . 

m,c 

. Uy{Q) . 

^xx ^xy 


UQ) ] 

. ByX Byy . 

m,c 

. ty{Q) . 


(3.15) 
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Taking each node in turn as the load point p and performing the integrations indicated 
in the above equation, a set of algebraic equations emerges as follows; 

[A]M = [5][t] (3.16) 

where the matrices [A] and [B] contains the integrals of the kernels T.j and Utj, 
respectively. Note that the parameter C,j{p) contributes only to the diagonal coefficients 
of the [A] matrix (when p is equal to Q). The formation of the matrices [A] and [B] can 
be seen as 


[A] 11 [A] 12 [A] 13 [A] 14 .... 


Ml 1 

[A] 21 [A] 22 [A] 23 [A] 24 •••• 


M 2 

[A]3i [A]32 [AJss [A]34 .... 


Ma 

[A] 41 [A] 42 [A] 43 [A] 44 .... 


M 4 

[B]n [B]u [B]i3 [5]i4 ■••• 


Wi ' 

[B]21 [B]22 \B]23 [-6] 24 •••■ 


W 2 

[B]3\ [B]z2 [.5]34 •••• 


[^]3 

[5]41 [B],2 [^]43 [BU •••• 


[tU 



m ' ■ * ' ^ 


where the submatrices of [A] and [B] are defined as follows: 


Mo = 


■^xx-^s^y 

• — 

BxxBxy 

AyxAyy 

J I-^JU 

ByxByy 


(3.18) 


Note that vectors [«] and [t] represent displacements and tractions in x and y direction 
as follows: 


[u], = 




(3.19) 


It is to be noted that the fundamental solution is ‘singular’. This singularity means 
that the kernels contains terms of the order ^ and fn(^) as 77 -+ 0. Since the kernels are 
dependent on the distance between p and Q, the three possibilities of the positions of p 
and Q can be examined as follows; 

p and Q in different elements: The standard Gaussian quadrature can be easily 
applied to the Uij and T,j kernels because they are not singular in this case: 


•/-l 3=1 


(3.20) 
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where G is the total number of Gaussian integration points (usually 4) and the 
Gaussian co-ordinates with an associated weighting function Wg. 

p and Q are in the same element but p Q: In this case, the and T,j kernels 
are singular but the shape function Nc{() in the vicinity of p is of the order r{p,Q). 
Therefore, the product of the kernels and the shape function is not singular and the 
integrals can be evaluted using the standard Gaussian quadrature. 

p and Q are in the same element and p = Q: In this case, the standard Gaussian 
quadrature can not be used because of the singularity of the kernels. Dealing with the 
C/,j kernels first, it is clear that as p coincides with Q, the singularity is of the form In 
(^) as p — » 0. This form of integral can be calculated by using the special logarithmic 
Gaussian quadrature and changing the limits of integration to become 0 to 1, as follows: 

/•i 1 

f{riVn{~)d7] = fiVgi)wgi (3.21) 

where Gi is the total number of logarithmic Gaussian integration points (usually 4) 
and rjgi is the Gaussian co-ordinate with an associated weighting function Wgi. A 
simple linear transformation can be used to accommodate the 0 to 1 integration limits 
as follows: 


1. If p is the first node of the element, rj = 0.5(1 + ^). 

2. Ifp is the second node of the element, the element is divided in to two 
sub-elements: rj = —( (for — 1 < ^ < 1). 

3. If p is the third node of the element, 7 / = 0.5(1 — (f). 

STEP 3: Application of the Boundary Conditions 

So far, calculation of all the coefficients of the [A] and the [B] matrices have been, 
discussed, but the problem is not unique until all the boundary conditions are applied. 
In a typical elastostatic problem, three types of boundary conditions are possible: 


1. Prescribed displacement. 

2. Prescribed traction (or stress). 

3. Linear relationship between traction and displacement. 

In a given problem, because there are 2N equations, it require 2N prescribed values. In 
other words, each node must have two of the four variables (ui, Uj,, fj., and ty) prescribed. 
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If a node has no prescribed value of any kind then it is assunaed to have = ty = ^ 
(i.e. stress free). 

To be able to use standard solvers, the [A] and [B\ matrices must be rearranged such 
that all the known variables are on the right-hand side and all the unknown variables 
on the left-hand side which will result in the following system of algebraic equations; 

[A*][;r] - [B*][y] (3.22) 

where [x] contains all unknown variables (whether displacement or traction), whereas 
[y] contains all known variables prescribed as boundary conditions (whether displace- 
ment or traction). The matrices [A*] and [B*\ are the modified forms of [A] and [B], 
respectively. Therefore, the final system of linear algebraic equations can be written as 
follows: 


[A*][a:] = [c 


where [c] now contains known coefficients. 


(3.23) 


To perform this arrangement of the [A] and [B] matrices, all the prescribed tractions 
are multipied by the relevant [B] coefficients and added together to form the relevant 
elements of vector [c]. If a node has a prescribed displacement then its [A] coefiicient is 
multiplied by the prescribed displacement value and moved to the right-hand side after 
reversing its sign. In its place, the [5] coefficient is moved to the left-hand side after 
reversing its sign. 

After solving the system of linear equations, the variables contained in the unknown 
vector [a:] of Eqn. (3.25) have to be sorted into displacements and tractions for computer 
output purposes. The above procedure can be easily applied to any size of mesh. 

Before solving the equations, it is worth noting that in practical elastostatic problems 
the displacement values are usually several orders of magnitude less than the traction 
magnitudes (because of the very high value of Young’s modulus, e.g. 200 xlO^ GPa 
for steel). This means that the typical value of the [A] coefl&cients are several orders 
of magnitude larger than the typical [B] coefficients. Since the [A] and [J5] coefficients 
are placed together in the [A*] matrix before solving the equations, they should be 
roughly of the same order of magnitude (otherwise, the solution may suffer from some 
inaccuracies). This can be easily achieved by multiplying the right-hand side of Eqn. 
(3.16) by a suitable scaling factor, as follows: 
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where the scaling foctor, S, is defined as follows; 


E 

S = (3.25) 

where E is Young’s modulus and Lmax is the maximum distance between any two 
nodes. Therefore, the traction values are effectively normalized, and after solving the 
equations can be simply multiplied by S to yield absolute (un-normalized) tractions. 

STEP 4: Solution of Algebraic Equation 

The solution matrix [A*] resulting after the application of the boundary conditions 
is not symmetric and fully populated with non-zero coefficients. Because of this, the 
choice of solution solver is clear: Gaussian eliminations or any other direct technique. It 
should be noted that all properly defined elastostatic problems must result in properly 
conditioned matrices (particularly if the scaling factor is used correctly) and they only 
become ill conditioned or singular if there is a mistake some where ( such as two nodes 
having the same co-ordinates ). It is a good practice to check for ill-conditioning when 
BE formulation is implemented in a computer program. 

STEP 5: Calculation of the Boundary Stresses 

There are two alternative methods of calculating the boundary stresses. In the first 
method, the BIEs for the stresses at p are used together with the appropriate third-order 
kernels. In the second method, the stresses are calculated directly from the computed 
displacements and tractions by simple differentiation using the shape functions. 

In Eqn. (3.7), it has been defined the unit tangential vector m, which has two 
components and rriy in the x and y directions, respectively. Now define the local 
directions 1 and 2 as tangential and normal directions, respectively. Therefore, the 
local tangential component of the displacement vector Ui can be written in terms of the 
cartesian displacements as follows: 


UiiO = + ^yiO'^y 

Using the shape functions from Eqn. (3.5) this becomes 


c=l C=1 




(3.26) 


(3.27) 
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To obtain the tangential strain en, the above equation is differentiated with respect to 
the tangential direction as follows: 


^u(0 = 



dN^iO 


(Uj;)c TTlj; -|- 



(3.28) 


Note that the Jacobian, J(^), is used in the above equation because of the relationship 
between the local co-ordinate ^ and the boundary path . 


The local components of the traction vector, ti and ^ 2 , can be defined as the tangential 
and normal tractions to the surface, respectively as shown in Figure 3.3. Note that 
the convention used here is that positive tangential traction should point to the left 
of the outward normal. If a is the angle between the normal traction and the global 
^-direction, then the local tractions can be written in terms of the cartesian global 
tractions as follows: 


ti = —txsina -f tycosa 

t2 = txcosa + tysina (3.29) 

To obtain the stresses in the local directions 1 and 2, the stress-strain are used as 
follows: 


<722 = h (3.30) 

cri2 = h 

Note that plane-strain conditions have been used in the above expressions for the 
stress. Provided that the ‘effective’ material properties of Eqn. (3.35) are used, the 
Eqns. (3.34) are also valid for plane-stress conditions. To transform local stresses into 
global stresses, the transformation matrix is to be used. The direction cosines of the 
local co-ordinates are easily obtained as 

E* — E\ u* = u] n* = {plane — strain) 

jg* = 7T~~! /^* = M (plane - stress) 

(l-f-i')® l + u 

Direction cosines of tangential direction 1 = (— sma); (cosa) 
Direction cosines of normal direction 2 = (cosa); (sina) 


(3.31) 

(3.32) 
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Figure 3.3 Local and Cartesian components of the traction vector. 


Therefore, using the transformation matrix ([4]), the global cartesin stresses are given 

by 

(Txx siri^a cos^a — 2sinacosa [ ctu 

(Tyy = cos^a sin^Q. 2sinacosa cr 22 (3.33) 

(^xy J L —sinacosa sinacosa (cos"^a — sin'^a) j [ . 

The angle can be witten in terms of the unit normal components, rix and Uy, as follows: 

a = (3.34) 

STEP 6: Calculation of Internal Variables 

After solving the equations, all of the values of displacements and tractions are de- 
termined at all the nodal points on the boundary. It is now a straightforward matter to 
determine the displacements at any interior point by substituting all of the boundary 
values in the BIE for displacements in Eqn. (3.1) and assuming the interior point to be 
the load point p. Similarly, the internal stresses can be determined by using the BIE 
for the stresses ([4]). 
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3.3 Sensitivity Analysis 


In this Section the sensitivity analysis has been carried out, in general, for any function 
with respect to the intrinsic shape design variables. The function may represent either an 
objective function or behavioural constraint of an optimization problem. For example, 
the objective function in case of shape optimization problems using BEM stress analysis 
could be the peak value of the Von-Mises equivalent stress at a particular node or it 
could be the weight of the component. The nature of the behavioural constraint will 
depend on the formulation of the optimization problem. Since the objective function 
or the constraint function is expressed in the cartesian domain, the sensitivity analysis 
requires the use of chain rule and change of variables concepts of partial differentiation. 

Consider a function g. It is assumed that the function g is explicitly related to the 
cartesian co-ordinates of the nodes, x„y,,i = 1,2, ,m. Each of these cartesian co- 

ordinates in turn related to the intrinsic parameters Kj,Sj,j = 0,1,2,- •• ,n. The rate 
of change of function g with respect to any intrinsic parameter can now be written as 

^9 — ^ I, ^ 1 

Oki ^^dxJ^dKi^ 

dg ^ A/ dg\ /dx,\ 

dsi ^^dxi^^dsi^ 

ingeneral, 

dg ^ ^/dg\^dxis 

dg _ ^■y/ dg \/ dxt \ 
dsn ^^dxj^dsj ^^^dyi'^dsn' 

The equations given above provide derivation of the function g can be estimated as 
follows. 


(3.35) 

(3.36) 

) 

) (3.37) 



(3.38) 


The foregoing procdure can be used to compute the sensitivity coefficients of any 
function g with respect to shape design variable. Also, Equation (3.42) can be used to 
compute the total change in a function g due to simultaneous changes in the SDVs. 


In the present work, the partial derivatives have been computed using the forward 
finite difference scheme. Some typical cases are given below. 
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^3 ^ g(a:i + A3 ;i,X2)- • • yi^y2,- ■ • ,ym) ~ gixi,X 2 ,- • • ,Xm, yi,y2,- ■ ■ ,ym) 

dxi Axi 

99_ ^ ,Xm, yi-h Ayi,y2, • • • ,ym) — g{xi,X2, • • • , , Vm) 

dyi Ayi 

(3.39) 

30,31,^2.’ ..Xri. 3o,3i,S2, ,an) 

A«1 

(^0 >^2 »• ap ,5^ 4* A^i ,32 (^0 t^il ,^2 , , a0iai,32» ‘M^n) 

9ai Aai 

(3.40) 

^Ul. 2:^ Vi(^Qi«1+A«i ,^2t jfiny S0,Si,S2, ^•iSn)-yi(K0,Ki,K2,^ ao,ai, 32 ....,an) 

dni A/ci 

9yi ^ ao,3i-f Aai ,32t"-ian)-Vl(Ko,/Ci,K2, -iKni ap ,3i ,32 i • lan) 

dsi Ail 

(3.41) 

It should be noted that appropriate values of Axi, Aj/i, A/Ci, and Asi should be 
chosen while computing the partial derivatives. In the present case a range of 0.1% to 
1.0% has been used to select a suitable value for the change of all these parameters. 
Computer programs have been developed in C on UNIX-platform to compute the sen- 
sitivity coefficients with respect to the intrinsic shape design variables. An illustrative 
example has been considered and the results obtained are given in the following Section. 

3.4 Numerical Examples 


A number of computer programs have been developed in both C and FORTRAN 77 
and integrated through system calls on SUN- workstation with UNIX-platform for shape 
synthesis, automatic mesh refinement of boundary elements, BEM analysis, NLP opti- 
mization and sensitivity analysis as explained in the foregoing Chapters. With the help 
of these computer programs the intermediate shape of an elastic ring under diametral 
compression problem has successfully been implemented. Once, the intermediate shape 
is achieved, a small perturbation on intermediate shape design variables ha.s been ap- 
plied and the resultant stress sensitivities are plotted a.s shown in Figure 3.7. 

Example; 


This example shows an elastic ring under diametral compression. Taking an advantage 
of symmetry, one quarter of the ring is modelled, as shown in Figure 3.4(a). Young’s 
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mudulus, Poisson’s ratio and the allowable stress are 2338 MPa, 0.38 and 18 MPa 
respectively. It is a case of plane-strain, and the object is to minimize the weight, 
while keeping the outer boundaries defining the exterior shape of the ring unchanged. 
For a shape optimization, T 3 and r4 are assumed to be fixed, r2 is allowed to move 
horizontally, and Pi is allowed to move freely. The boundary Pi is being varied by 
selecting a feasible shape model such that stress constraints are satisfied. 

Figure 3.4(b) shows quadratic boundary element model with 100 nodes and 50 ele- 
ments. The cartesian and intrinsic profiles for an initial geometry of the ring as shown 
in Figure 3.5. A 2-LINCEs model (R-R model) has been chosen in order to define each 
boundary Fi, r2, Fs, and F4. The intrinsic variables of the boundary Fi (i.e. curve 
A-B) corresponding to the intial geometry are given by 

«o = -26.315789; sq = 0.0; /Cj = -26.315789; si = 0.0298451; K2 = -26.315789; 
52 = 0.0596902 

The Figure 3.6 shows the cartesian and intrinsic profiles for an intermediate shape of 
an elastic ring. The intrinsic shape design variables of the boundary Fi corresponding 
to an intermediate shape are given below: 

Ko = -26.315789; sq = 0.0; = -19.614253; Si = 0.029845; = -35.08690; 

52 = 0.062217 

In order to calculate the stress sensitivity with respect to Ki, a set of 20 different 
magnitudes of perturbed values (i.e. A«i) have been applied to the intermediate value 
of Ki. The Figure 3.7 shows the stress sensitivity Vs change in curvature according to 
the stress sensitivity procedure described in Section 3.3. This stress sensitivity study 
gives an information that a small change in intrinsic variables (i.e. curvature or arc 
length) caused a large variation in stress magnitude. 




Figure 3.4 (a) Definition of An Elastic Ring Problem. 

(b) BEM Mesh for the Elastic Ring Problem 
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Stress Sensitivity 
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Chapter 4 


SHAPE OPTIMIZATION 
METHODOLOGY 


4.1 General Approach 


The first and foremost step towards shape optimization based on a CAD approach 
is the creation of a multi-window interactive graphics environment in order to facilitate 
clear visualization of the design history of the optimal design process. To fulfil this 
requirement, a computer program has been developed in C language and implemented 
on a SUN-workstation with UNIX-platform based on SUN-view graphics. The graphics 
environment is comprised of multi-windows; one for the cartesian space, one for the 
intrinsic space and one for the selection menu and the last one for the graphics manipu- 
lation during the design process. Figure 4.1 shows such a graphics environment. Once, 
such a graphics environment is ready, then the shape optimization methodology can be 
described as follows. 

The problem of shape optimization can be described as a problem of synthesizing 
the geometry of the surface boundary of a component cross-section that will maximize 
or minimize an objective function as well as satisfy a set of behavioural and side con- 
straints. In general, the equation of the objective function and constraints require the 
geometry of a component surface boundary in terms of cartesian co-ordinates. How- 
ever, as proposed in this research, the shape definition is described most suitably in 
terms of intrinsic geometry of the component. The structural and mechanical compo- 
nents considered for this research use BEM analysis for 2-D plane-stress, plane-strain 
and axi-symmetric problems in order to satisfy the stress constraints during optimiza- 
tion cycle. In order to avoid the intervention of the designer during shape optimization 
cycle, an automatic mesh refinement of boundary elements has been introduced. The 
automatic mesh refinement module acts as a mediator between the synthesizer and the 
BEM analyzer. 
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Figure 4.1 Multi-Window Graphics Environment Designed 
on SUN“Workstation via SUN-view Graphics. 
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Consider a given component F (surface boundary) = Fi U r 2 U Pa U r 4 as shown in 
Figure 4 2. It is required to find a set of optimal shape curves between successive pairs of 
points such as Fi, Fa, F 3 ,andF 4 . Let ■ ■ • ,*n) be a set of n shape design variables 

(SDVs). These are considered to be intrinsic geometry parameters. Let (ci,C 2 , . . . ,Cm) 
be a set of m variables representing the cartesian geometry of the component. Let M 
be a mapping relation that enables us to find the set {ci, c^, . . . ,Cm) from a given values 
of ■ • • ,*n)- 

(<^1) • • • , <^m) — [-^](h) *2) • ■ • , in) (4-1) 

Let W be the objective function and g/s be the set of constraints. The optimization 
poblem formulation can be stated as follows. 

Minimize 


fF(ii,i2, ,in) (4.2) 

subject to 

t7j(^l)^2, )*n)) (^•^) 

where j = 1 , 2 , ...., fc. 


At this point, it is necessary to point out that Ci,C 2 , ,c„i are not the design variables 

for the present problem. Any change of the shape is effected by changing the vector 
(il ) i2j ) in)- 

The design poblem, in turn, now is to select a shape model for ri,r 2 ,r 3 , and F 4 etc. 
individually or collectively and their shape design variables such that W is minimized. 
The 2 -LINCEs model (R-R model) and 3 -LINCEs model (R-R-R model) have been con- 
sidered in this research work, but the same procedure can be extended to n-LINCEs 
model depending upon the problem situation. In general, the R-R-R model ( 3 -LINCEs 
model) has been found to be suitable having a minimum number of shape design vari- 
ables. The solution procedures via zero-order and first-order NLP optimization methods 
outlined in Sections 4.2 and 4.3. The sensitivity aspects have been discussed in Sec- 
tion 4.4. Before outlining the zero-order method, it is necessary to understand how the 
frontal-area of the surface boundary can be evaluated, if the geometry is available in B?- 
Euclidean space. The following equation has been used for evaluating the frontal-area 
of the closed boundary in a 2 -Dimensional space. The weight of a component can be 
evaluated once the cross-sectional area is available. 

Area = f. (4.4) 

t =:0 ^ 

where the vertices of the closed polygon are marked as 0 , 1 , 2 , • • • , (n — 1). The starting 
point is also designated as n. 
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Figure 4.2 General Definition of the Shape Optimization 
Problem Using Intrinsic Geometry and BEM. 
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4.2 Zero-order Method (Exhaustive Search Method) 


It is convenient to categorize algorithms according to the type of information that 
must be provided in searching for the minimum of the function. 

The simplest approach to minimizing IT(X) is to select either randomly or in some 
logical manner, a large number of candidate X vectors and evaluate the objective for 
each of them. The X corresponding to the minimum W{X) obtained from this set is 
called the optimum, X*. Obviously, if a precise solution to the problem is to be found, a 
great many X vectors may have to be considered. Methods such as this, which require 
only function values in searching for the optimum, are referred to as zero-order methods. 

These methods are usually reliable and easy to program, often can deal effectively 
with nonconvex and discontinuous functions, and in many cases can work with discrete 
values of design variables. The price paid for this generality is that these methods often 
require thousands of function evaluations to achieve the optimum, even for the simplest 
problems. Therefore, these methods are considerd most useful for the problems in which 
the function evaluation is not computationally expensive. Used in appropriate fashion, 
zero-order methods are a useful and powerful design tool for many common applications. 

Exhaustive Search method, in particular can be used to solve problems where the 
interval in which the optimum known to lie is finite. Let x, and xj denote, the starting 
and final points of the interval of uncertainity. The exhaustive search method consists 
of evaluating the objective function at predetermined number of equally spaced points 
in the interval (x^, x/), and reducing the interval of uncertainity, assuming that the 
objective function is a unimodal function. In general, if the function is evaluated at n 
equally spaced points in the original interval of length Lq = Xf — x, and if the optimum 
value of the function (among the n function values) turns out to be at the point xy, then 
the final interval of uncertainity is given by ([41]) 

Ln = Xj+t - Xj_i = (— -— )Lo (4.5) 

71 'Y JL 

It is to be noted that the objective function may be defined in terms of intrinsic 
variables explicitly or in terms of cartesian co-ordinates of some specific interpolating 
points. For the applications dealt in the present work the objective function involves 
in many cases the computation of the frontal-aera. The frontal-area in turn, can be 
calculated using the formula given under Equation (4.4) earlier. Moreover, the geometry 
or the curve is synthesized using a 2-LINCEs or 3-LINCEs model. 

The shape optimization methodology via a zero-order method can be described as 
follows. 
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Minimize 


subject to 


$ = $ 


Cl,C2, 


i(*l) *2 



(4.6) 


9j 


Cl,C2, • • 


i(*i,*2r--,*n)j - gall < 0 


(4.7) 


where j = 1,2, 


The description of the terms involved in the above two equations include: $ defines 
the frontal-area of cross-section of the component, the outer boundary of the component, 
and gj defines the peak value of Von-Mises stress corresponding to the jth iteration, and 
gall defines the maximum allowable stress of the component. The solution procedure 
can be outlined as follows. 


STEP 1. Define the surface boundary (P =ri Ur 2 UPa Ur 4 etc.) in terms of intrinsic 
parameters by selecting a shape model depending on the physical nature of the problem. 
In general, the 3 -LINCEs model (R-R-R model) has been found to be suitable having a 
minimized number of shape design variables (SDVs). 

STEP 2 . Once a particular model is chosen, the SDVs are selected depending on the 
number of free parameters for that particular model. It has been found that s serve 
as a better choice for shape design variables than s,' s. 

STEP 3. Define cartesian-intrinsic mapping for a feasible set of SDVs, find cartesian 
co-ordinates using equations (2.31) through (2.33). 

STEP 4. Once the cartesian co-ordinates are available, it is possible to analyze stress 
constraints (g/s - gaii < 0) for maximum Equivalent stress (Von-Mises stress), maximum 
deflection etc. using Boundary Element Method. If any constraint is violated then the 
design is discarded. The process now resumes at step 3 with a new set of shape design 
variables. This step also involves the automatic mesh refinement of boundary elements 
through an automatic mesh refinement module. 

STEP 5. Check how stress and deflection situation changes by altering Pi, and hence 
P 2 ,P 3 , and P 4 individually or simultaneously. 

STEP 6 . Evaluate objective function (i.e. frontal-area $ in this case) and optimal 
shapes using values of shape design variables and the cartesian co-ordinates that satisfy 
all the constraints. If the value of ^ is minimum then the process terminates. Otherwise 
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an improved feasible design is required. This will be accomplished by an exhaustive 
search NLP (Non-Linear Programming) technique. 

The exhaustive search NLP technique which has been used in this method can be 
summarized as follows. Initially the interval of uncertainity of free SDVs is assumed 
to lie in between K^ax and /Cmtn corresponding to inner and outer boundaries of the 
frontal-area of component being minimized respectively. The interval of uncertainity 
will be reduced with the assumption of unimodal function and the accuracy with which 
the optimum value is required. This decides the number of iterations. Once the number 
of intervals is decided based on the above procedure, it is possible to select a set of 
values of shape design variables - some s, s and s. The values of remaining /c/s and 
s.'s are obtained from the geometric constraint equations. The SDVs thus obtained will 
be used in establishing the shape model and hence the objective function $. 

The shapes which satisfy all the stress constraints can only undergo the frontal-area 
evaluation whereas the shapes which are not able to fulfil any of the stress constraints will 
be discarded. The value of $ which will be obtained in every iteration of the optimization 
cycle is compared to the existing or preceding one for the condition of minimization. 
Thus the shape corresponding to an optimal frontal-area can be obtained. 

4.3 First-order Method 


Let VP(X) be the objective function, then a more difficult but usually more efficient 
approach to the minimization problem is to use gradient information in seeking the 
optimum. For example, it can be said that calculate the gradient of W (X) and then 
search in the negative, or — VW(X), direction. Now because this is the direction of 
steepest descent at X, then it may be expected that choosing a new X vector using this 
information will lead to move to a solution of the problem rapidly. That is, a new X 
vector will be chosen by 


X^ = X° - aVlF(X°) (4.8) 

If several values are tried for the scalar parameter a, and for each resulting X vector 
evaluate VF’(X) = kF(Q:), then some value a* will provide a minimum 14^(X) for this 
search direction. By using the gradient of 1F(X), it will be possible to limit our search 
to a specific direction, rather than randomly searching the entire space. Of course, 
because VF(X) may be highly non-linear, then it will be necessary to calculate a new 
gradient at X and repeat the process. Methods such as this which use gradient, or first 
derivative, information are called first-order methods. It is to be noted that the objective 
function and geometric and stress or behavioural constraints must be expressed in terms 
of intrinsic shape design variables either explicitly or implicitly. The formulation has 
been done for the problems which are characterized 2-D plane-stress or plane-strain case. 
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A 2-LINCEs model (R-R model) and a 3-LINCEs model (R-R-R model) have been used 
for modeling the surface boundary of the component. 

The problem formulation (for 3-LINCEs model) can be stated as follows. 

Minimize 


^ ^ *^2) ‘ ’ 5 ^2? ' ' ' 5 *n)| 

subject to 

(i) Geometric Equivality constraints (h, s) include : 


(4.9) 


X3 - Xq 


I 

I 


COS 


S2 


COS 


Si 

^3 


J32 


COS 


^)(^ " soa) + kqct + Cl 

{— ^)(^ - ^1<^) + + <^2 

L 32 ^ Z 

(— ^)(^ - 520^) + «2<7' + C3 

1^3 — S2 2 


da 

dor 

dcr = 0 


(4.10) 


ys - yo 


rsin\{ 

JsQ 

r ■ 

I sin 
J 3l 


' ^ ^)(’T — •^oor) + /CoCT + Cijrfcr 

3i — Sq Z J 

- Sicr) -h Ki<7 + C'^da 

,52 3i Z J 


(4.11) 


— / 3in[(”^^ — 52^t) + + C'a — 0 

J32 L ’Ss 32 " 

^n-V’O- ~- -^(3i - So) - -5l) - (^3 “ • 82 ) = 0 (4.12) 

(ii) Inequivality stress constraint [(gj — Qaii) s] and other geometric constraints include : 


{Tjnax Tall) ^ 0 


^max 1 0 


(4.13) 


or 
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^3 ^ ^mtn 


or -X 3 - < 0 


(4.14) 


^3 ^ ^max 


or X3 - Xmax < 0 

(4.14 a) 

(iii) Constant values include : 

If t^tart and -i/jend = Specified tangent conditions at the start and end of boundary 
curve, 


then ?/>o = V’aiart 
^3 = i^end 

xo = 0 (4.15) 

J/3 = 0 

and So = 0 

The terms involved in Eqns. (4.9) through (4.15) have got their own meaning as 
described in Chapter 2 and Chapter 3. T^ax corresponds to the peak value of Von- 
Mises equivalent stress at the boundary and Tail defines maximum allowable stress of 
the component. 

In order to convert the above constrained non-linear optimization into unconstrained 
optimization problem, the following SUMT (Sequential Unconstrained Minimization 
Technique) has been used. 

If I4^(X) is an objective function, then 

minimize 


PV'(X) 


(4.16) 


subject to 
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5fj(X) < 0 w/iere j = 1,2, • • • ,m (4-17) 

and /i,(X) = 0 where i = 1,2^- ,p. (4-18) 

Using interior penalty function method, the above constrained optimization problem 
can be converted into unconstrained minimization problem by adding suitable penalty 
terms. The pseudo-objective function after adding suitable penalty terms using interior 
penalty function can be written as follows ([41]). 

Pseudo-objective function include ; 


= $(X,rfc) 

mi IP 

where is called the penalty parameter. 

The solution procedure for the proposed first-order shape optimization methodology 
can be outline as follows. 

STEP 1. Define the surface boundary (P = Pi U r2 U Pa U P4 etc.) of the initial 
geometry in terms of intrinsic parameters by selecting a shape model depending on the 
physical nature of the problem. In general, the R-R-R model has been found to be more 
suitable having less number of shape design variables. 

STEP 2. Once a particular shape model is chosen, the SDVs are selected depending 
on the initial geometry of the problem. 

STEPS. Define intrinsic-cartesian mapping for a feasible set of SDVs, find the carte- 
sian co-ordinates using shape synthesis procedure explained in Chapter 2. 

STEP 4. An automatic mesh generation of boundary elements will be obtained using 
an automatic mesh refinement module, once the cartesian co-ordinates of geometry, 
displacement and stress boundary conditions and material properties are available. 

STEP 5. Evaluate the maximum Von-Mises equivalent stress (xmox) at the surface 
boundary using Boundary Element Method. 

STEP 6. Since the maximum allowable stress (tj;;) of the component is given, then the 
first-order constrained optimization problem can be formulated using Equations (4.&) 
through (4.15) and then normalizing the necessary constraints. 

STEP 7. Once the Non-linear constrained optimization problem is formulated as 
stated in STEP 6, then the non-linear constrained optimization problem will be con- 
verted in to a non-linear unconstrained optimization using interior penalty function 
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method of Secjuential Unconstrained Minimization Technicjue (SUMT) according to the 
Equation 4.19. Thus, the non-linear constrained problem becomes non-linear uncon- 
strained optimization problem with suitable penalty parameters added to the objective 
function, called pseudo-objective function. 

STEP 8. Once the non-linear unconstrained optimization problem is formulated, 
then the optimization solution can be obtained using Davidon-Fletcher-Powell (DFP) 
method of unconstrained minimization technique. The gradient information about the 
objective function required to be used may utilize concepts of chain rule and change 
of variables for partial differentiation depending upon the objective function is defined 
either explicitly or implicitly in terms of the SDVs. 

STEP 9. The optimal shape of the component can be visualized using the graphics 
environment created on SUN-workstation as mentioned at the outset of this Chapter 
with the help of cartesian geometry produced in STEP 8. 

STEP 10. Using the cartesian geometry corresponds to an optimal shape evaluate 
the objective function (i.e. the frontal-area of the cross-section) of the problem. Thus, 
the shape optimization cycle will be accomplished. 

4,4 Issues of Sensitivity Analysis 

The sensitivity analysis is an approach to indicate how a constraint behaves owing to 
small perturbations in the values of the design varibles. In an intrinsic geometry based 
shape optimization, the effects of small changes in intrinsic shape design variables, Ki,si 
etc, on the behaviour of constraints will be studied in this section. For example Figure 
4.3(a) and 4.3(b) shows cartesian and intrinsic geometries respectively. The sensitivity 
analysis concept based on intrinsic geometry can be explained with the help of these 
Figures. 

If there is a change in the value of, say Ki to /Cj -f- A/ci then the effect on the x-y curve 
(curve A-B) can be studied. 

Suppose there is function such as an area (A) which depends on the x-y (and hence 
on K-s shape) (of the) curve, then the change 

A A in the area (A) due to a change in the variable «i as A«i. 

If Ki is a shape design variable and area (A) is a constraint then the ratio (^;:;]|') is 
the sensitivity of the area with respect to a shape design variable; in this case ki. 

If it is case for /ci then the same logic could be used to evaluate for any other SDVs such 
as «o, « 2 , ...., «„ or so,si, ...,, The solution procedure for sensitivity computations can 
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be described using chain rule and change of variables concept of partial differentiation 
and finite difference approximations. 

If $ = ^{kq, So, «i, Si, • • • , Kn, Sn) = $(I) defines the frontal-area of the surface bound- 
ary, then the rate of change of # w.r.t can be defined as follows. 

Where 


1 = 


Ko 

■30 

K2 

S2 


> 


«n 
Sn } 


then 


5$(/) _ $(I(/ci -f Aki)) - $(I) 
dni 


(4.20) 


where A/ci denotes a small perturbation in ki. Similarly, for other SDVs ko,so,K 2,32 
etc. can be defined. 


Suppose $ = $(ci,C 2 , = ^(^(I)) define the frontal-area of the 

surface boundary, then the sensitivity computations will be carried out as follows: 
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C = 


Xi 

yi 

X2 

2/2 


Vm 


Where (ci, C 2 , c^) defines cartesian variables such as (a:i,i[/i),(x2, J/ 2 ), ■•■,(a:m,ym) 
and ( 21 , 22 , defines intrinsic variables such as kq, ^o, «!, si, •••■,/««, 

Then 




(4.21) 

but 

yi,y2," yi,y2, -.ym) 

dxi 

(4.22) 


^(a7i,ar2r-,a;m, I/l-f Ai/i Jfl 


dyi 

Ai/i 

(4.23) 

dxi ^ 



dni 

Akj 

(4.24) 

and 

dy\ S/l («Oi«l 4* A/Ci ,«2, — M^n) 

dK\ Aki 

(4.25) 


Similarly, the same concept could be used for the other shape design, variables. 


If more than one shape design variable is perturbed simultaneously, for example k\ and 
Si then the net effect on the frontal-area (i.e. objective function) would be determined 
as follows. 



(4.26) 

The value of gives the total change of $ w.r.t the simultaneous variation of 
and 3 i. The same concept could be extended to n-shape design variables. This sensi^ 
tivity procedure has been implemented on a numerical example of an elastic ring under 
diametral compression in Chapter 5* 
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4.5 Comparison and Efficiency of Proposed Method 


As mentioned at the outset of this Chapter, the interactive graphics capabilities in- 
tegrated with the proposed shape optimization cycle is prudent to allow designer in- 
teraction in the Computer Aided Design (CAD) process. Such a dialogue can be very 
beneficial, saving computer and human resources. 

In general purpose design optimization methodology needs the following information 
about the problem to be solved : (i) input data such as number of design variables, 
number of constraints, etc. (ii) the cost and constraint functions (iii) gradient of cost 
and constraint functions. If the gradients are not available system should automatically 
approximate them by a finite difference method. If there is a mistake in the input data 
or problem definition, error will occur in the problem-solving procedure. The optimiza- 
tion system should take care of such mistakes as much as possible. In contrast, the 
interactive graphics capabilities developed on SUN-workstation and HP 9000/800-series 
workstation wth UNIX-platform together with intrinsic geometry based shape synthe- 
sis, automatic boundary element mesh generation and refinement, BEM analysis, and 
numerical optimization programming modules makes the shape optimization cycle more 
attractive. The introduction of multi-window concept makes the intrinsic shape defini- 
tion more versatile. This is due to the fact that the effect of changes in the intrinsic 
space can be seen immediately in the cartesian space, thus it gives real feeling of the 
shape to the designer. In addition, the interactive capabilities include : Design history of 
cost function, constraint functions, design variables, maximum constraint violation, con- 
vergence parameter would be monitored. Since, the interactive environment developed 
in this approach, especially graphical display can be of great help in certain decision 
making during design process. 

Since the computational cost of shape optimization is very sensitive to the number 
of shape design variables needed to define the boundary geometry, the number of SDVs 
can be kept to a minimum using intrinsic shape synthesis approach compared to Bezier 
or B-spline and Hermitian cubic spline functions. Hence, it can be said that intrinsic 
form provides flexible and compact geometric capability to the designer. 

The strategy that has been adopted in this research involves mapping the geometry 
from intrinsic to cartesian domain and then interpolating to quadratic boundary ele- 
ments for BEM stress analysis purpose. There is no direct relation between cartesian 
co-ordintes obtained after intrinsic-cartesian mapping and the boundary element nodes, 
although in some cases they may happen to coincide. 

The number of iterations required via zero-oder method in the proposed shape opti- 
mization cycle is slightly more than the number of iterations required in the reviewed 
methods of shape optimization using B-spline or Hermitian geometric representation 
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and BEM analysis. Wlieieas, the first-order proposed method requires less number 
of iterations when compared to the methods based on B-spline or Hermitian geomet- 
ric representation and BEM analysis. The present methodology tackles the linear and 
non-linear objective functions and side and behavioural constraints. 

The total cmpputational (CPU) time required for the proposed shape optimization 
cycle vary from six to seven hours on SUN-workstation depending upon the nature of 
the problem and the initial geometry. 

The proposed shape optimization methodology can be applied to a broad range of 
applications in order to show its effectivenss over the parametric approach. The com- 
puter programs developed based on the above mentioned shape optimization methodol- 
ogy include: Non-linear Equation solver, Numerical Integration module, An Automatic 
Boundary Element Mesh generation and refinement module, Boundary Element Analysis 
for Design (BEAD), and Zero-order (or First-order) NLP optimizer makes the software 
elegant for 2-D elastostatic problems characterized by plane-stress, plane-strain, and axi- 
symmetric cases. It deserves particular acclaim in its ability to get results in low number 
of cycles, few analyses of small dimensionality, and proper use of reliable optimization 
algorithms, capable of making best use of approximation concepts in minimizing the per 
step optimization problem size. 

Aside from the above built in efficiency properties, it is however, to be noted, that 
the proposed separation of shape synthesis module from the analysis model, allows the 
designer to make a judicious choice in defining the shape via simple intrinsic forms and 
low number of design variables, while at the same time using a highly refined analysis 
model for sufficient accuracy of the performance of the end results. Therefore, the 
matter of overall efficiency of the performance of the proposed CAD approach, remains 
to a great extent in the hands of the designer. 



Chapter 5 


CASE STUDIES VIA 
ZERO-ORDER METHOD 


5.1 Computer Software Development 


The intrinsic definition of a curve or a surface seems to suffer from a handicap that 
one doesn’t get the feel of the geometry explicitly. In order to eliminate this drawback of 
intrinsic geometry a multi-window platform has been developed using SUN-view graph- 
ics on SUN-workstation and STARE ASE graphics on HP 9000/800-series workstation. 
The programs have been written in C-language. In additon to the graphics environment 
created as said above, number of programs for non-linear equation solver, numerical 
integration, shape synthesizer, automatic boundary element mesh generation and re- 
finement, Boundary Element Analysis for Design (BEAD), and NLP optimizer have 
been developed in C and FORTRAN 77 languages. The linking between C and FOR- 
TRAN 77 has been carried out via system calls in Unix-environment. The illustrative 
examples described in this Chapter are modelled using either 2-LINCEs model (R-R 
model) or 3-LINCEs model (R-R-R model) and are characterized by 2-D elastostatic 
plane-stress or plane-strain case. Shape optimization is carried out for minimum weight 
design problems. Grapher and ANSYS/FEM packages have been used for plotting 
stress versus node number graphs and colour stress plots respectively. Also, number of 
C-programs have been developed based on sensitivity procedures explained in Chapter 
3 and Chapter 4. 

5.2 Illustrative Examples 


Based on the proposed shape optimization methodology via zero-order method (ex- 
haustive search method) practical engineering problems include: (i) A thick cylinder 
subjected to a constant internal pressure, (ii) an elastic ring under diametral compres- 
sion, (iii) a torque arm subjected to an axial and transverse loading, (iv) a fillet under 
an axial load, and (v) a ladle hook subjected to a tensile load have successfully been 
implemented. The sensitivity analysis has been carried out for an elastic ring problem. 
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$1 - Area bounded by the outer curve E-F (for a quadrant section) = f (radius of the 
outer boundary)^ 

Tmax - Peak value of Von-Miscs stress at the surface boundary. 

Tall - Allowable stress or Yield stress of the component. 

Optimal shape has been obtained using shape optimization algorithm explained in 
Section 4.2. From Figure 5.3(b) the shape design variables related to initial shape 
including: kq = —13.3333; so = 0.0; ki = -13.3333; sj = 0.0589048; K 2 = -13.3333; 
and 32 = 0.1178097. The Figure 5.4(b) shows optimal shape related to SDVs including: 
Ko = —13.3333; sq = 0.0; ki = —10.405382; sj = 0.04; K 2 = —16.134596; and S 2 = 
0.122595. The frontal-area (i.e. objective function) has been reduced 10.34% from the 
uniform thick cylinder, with the inner radius the same as the smallest in optimum. 

5.2.2 An Elastic Ring Problem 


This example shows an elastic ring under diametral compression. Taking the ad- 
vantage of symmetry, one quarter of the ring is modelled, as shown in Figure 5.9(a). 
Young’s modulus, Poisson’s ratio and the allowable stress are 2338 MPa, 0.38, and 18 
MPa respectively. The boundary Fi (i.e. curve A-B) is varied in the search direction 
by selecting a suitable shape model such that stress constraints are satisfied. The point 
B or C is assumed to be free to move horizontally in this problem, as shown in Figure 
5.9. The other boundaries Fa (i.e. curve E-F) and F^ (i.e curve F-G) are assumed to 
be fixed. The stress analysis is done for the plane-strain condition. Figure 5.9(b) shows 
quadratic boundary element model with 100 nodes and 50 elements. 

The shape model chosen in order to vary the curve A-B to be of 2-LINCEs model 
(R-R model) with Ki cind si as free shape design variables. The formulation of the above 
problem can be stated as follows. 

Minimize 


^ $ 


[/'" 

Lt/.to 


“ -So) +• «o} 


ds -1- 


Jsi 




ds (5.3) 


subject to 


'^max '^all ^ 9 


(5.4) 


where 

^ - Defines the frontal-area of the surface boundary of a component (for quadrant 
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Figure 5.3 (a) Cartesian Profile for the Initial Shape of 
a Thick Cylinder. 

(b) Intrinsic Profile for the Initial Shape of 
a Thick Cylinder. 
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ure 5.5 (a) Cartesian Profiles of Initial and Optimal 
Shapes of a Thick Cylinder. , 

(b) Intrinsic Profiles of Initial and p 
Shapes of a Thick Cylinder. 
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Figure 5.6 Cartesian and Intrinsic Profiles of Initial 
and Optimal Shapes of a Thick Cylinder 
on a Multi-Window SUN-view Graphics 
Environment. 




Equivalent Stress (in Pa) 



Figure 5.7 Equivalent Stress Vs Node Number along the 
Boundary of a Thick Cylinder. 
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$1 - Area bounded by the outer curve E-F (for a quadrant section) = |(radius of the 
outer boundary)^. 

Tmoi ■ Peak value of Von-Mises stress at the surface boundary. 

Tall - Allowable stress or yield stress of the component. 

Optimal shape has been obtained using shape optimization algorithm explained in 
Section 4.2. From Figure 5.10(b) the shape design variables related to initial shape 
including; kq = -26.315789; sq = 0.0; «i = -26.315789; 5 i = 0.0298451; = 

—26.315789; and S 2 — 0.0596902. And the Figure 5.11(b) shows SDVs related to an 
optimal shape include : kq = -20.315789; 5o = 0.0; ki = -16.797924; sj = 0.0120; 
K 2 — —33.063226; and S 2 = 0.078075. The frontal-area (i.e. objective function) has 
been reduced 17.9% from the uniform ring section, with the inner radius the same as 
the smallest in optimum. 

Figure 5.10(a) and (b) shows cartesian and intrinsic profiles for an initial shape of an 
elastic ring. Cartesian and intrinsic profiles of an optimal shape are shown in Figure 
5.11(a) and (b). The cartesian and intrinsic profiles of an initial and optimal shapes 
together are included in the Figure 5.12(a) and (b) respectively. 

Figure 5.13, show design history of an elastic ring in SUN- view graphics environment. 
Cartesian and intrinsic geometries of an elastic ring corresponding to initial and optimal 
shapes in a multi-window SUN-view graphics environment has been shown in Figure 
5.14. 

Von-Mises Equivalent stress versus Node number graphic plots of initial and optimal 
shapes are shown in Figure 5.15. The colour stress plots belong to initial and optimal 
shapes of an elastic ring problem have been shown in Figure 5.16. 

5.2.3 A Torque Arm Problem 

The torque arm in Figure 5.17(a) has been optimized for minimum weight (i.e. min- 
imization of frontal-area) according to the shape optimization procedure described in 
Section 4.2. The hole shown at the bigger end is where the device is fastened to the 
support taken with zero displacements, while the applied load is extended at the hole 
in the smaller end where force boundary conditions are given. The magnitude of the 
applied loads, allowable stresses and material properties are also given in Figure 5.17(a). 

The objective is to minimize the weight of the arm by keeping the thickness of the 
material constant and finding the most suitable shape offering proper distribution of 
material, for the most efficient load carrying capacity. 

Clearly, the problem is a plane-stress case and it is modelled with quadratic boundary 
elements. The whole boundary is divided into 80 elements and 160 nodes. Figure 5.17(b) 
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Figure 5 9 (a) Definition of An Elastic Ring Problem 
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Figure 5.10 (a) Cartesian Profile for the Initial Shape 

of An Elastic Ring. 

(b) Intrinsic Profile for the Initial Shape 
of An Elastic Ring. 
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Figure 5.12 (a) Cartesian Profiles for the Initial and 

Optimal Shapes of An Elastic Ring, 
(b) Intrinsic Profiles for the Initial and 
Optimal Shapes of An Elastic Ring. 
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Figure 5.13 Design History of An Elastic Ring on a 
Multi-Window SUN-view Graphics 
Environment. 




Chapter 5. Case Studies Via Zero-order Method 






Equivalent Stress 


Chapter 5. Case Studies Via Zero-order MethnH 







(Chapter 5. Case Studies Via Zero-order Method 


108 




Chapter 5. Case Studies Via Zero-order Method 


109 


shows the boundary element mesh. 

Allowing for the shape variations in the mid-section of the arm, located between the 
holes (i.e curve C-D and curve G-H). These two curves are modelled using 3-LINCEs 
model (R-R-R model) for each of the curves with suitable shape design variables. A limit 
on the maximum Von-Mises stress in an element is imposed. It is also to be noted that 
the modeling is to be done by keeping the holes and the semi-circular regions around the 
holes unchanged. The intrinsic SDVs related to an initial shape can be seen from the 
Figure 5.18(b) and the intrinsic SDVs related to an optimal shape are shown in Figure 
5.19(b). 

Starting from the shape given 5.18(a), and using the criterion for design as shown 
below. 

Minimize 


$ $(ci , C 2 , ■ ■ • , C,7j (il , *2, ■ ■ ' , *n)) 


(5,5) 


subject to 


'^max Tall ^ 0 (^■®) 

the design undergoes the process of changes according to the shape optimization method- 
ology described in Section 4.2 and finally gives an optimal shape of the torque arm as 
shown in Figure 5.19(a). Figure 5.20(a) shows the cartesian profiles corresponding to 
initial and optimal shapes of the arm. The intrinsic profiles corresponding to initial and 
optimal shape of the arm are shown in Figure 5.20(b). The total of 45.6% reduction in 
frontal-area has been achieved. 

In Figure 5.21 the cartesian and intrinsic profiles for initial and optimal shapes of the 
torque arm in a Multi-Window SUN-view graphics environment has been shown. The 
colour stress plots of initial and optimal shapes are shown in Figure 5.23. 

5.2.4 A Fillet Problem 


It is required here to find the best shape of the fillet that makes a transition between 
two thicknesses of plate material connected in a butt joint. Figure 5.23(a) shows the 
problem with the tensile applied load and the material properties. A limit on the 
maximum Von-Mises stress in an element is imposed (Tmax)- The objective function 
would be the area of material enclosed in the triangle ABC to be minimized. The BEM 
model is setup with 160 nodes, with 80 quadratic boundary elements, as shown in Figure 
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Figure 5.17 (a) Definition of a Torque Arm Design 

Problem with Displacement and 
Stress Boundary Conditions, 

(b) Boundary Element Mesh for a Torque 

Arm. 








Figure 5.18 (a) Cartesian Profile for the Initial Shape 

of a Torque Arm. 

(b) Intrinsic Profile for the Optimal Shape 
of a Torque Arm. 
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Figure 5.20 (a) Cartesian Profiles for the Initial and 

Optimal Shapes of a Torque Arm. 

(b) Intrinsic Profiles for the Initial and 
Optimal Shapes of a Torque Arm. 
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Figure 5.21 Cartesian and Intrinsic Profiles for the 
Initial and Optimal Shapes of a Torque 
Arm on a Multi-Window SUN-view Graphics 
Environment. 
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5.23(b). The applied load consists of 1777 N uniform on the cross-sectional net area. 
To define the shape of the transition section between fixed points A and B, a 3-LINCEs 
model (R-R-R model) is being used. 

Formulation of the problem can be described as shown below. 

Minimize 


subject to 


$ ^(ci, C2, • * • j j ^2) ‘ * 5 ^n)) 


(5.7) 


'^max 

'^all 


-1 <0 


(5.8) 


where 

# - Frontal-area of the total cross-section of the plate-fillet joint assembly. 
Tmax - Maximum Von-Mises stress at the surface boundary. 

Tall - Allowable stress of the fillet joint. 


Figure 5.24(a) shows cartesian profile and Figure 5.24(b) shows intrinsic profile for an 
initial shape of a fillet. Cartesian and intrinsic profiles for an optimal shape are shown 
in Figure 5.25. Combined cartesian and intrinsic profiles for an initial and optimal 
shapes are shown in Figure 5.26. A multi-window SUN-view graphics environment with 
cartesian and intrinsic profiles of a fillet is shown in Figure 5.27. Colour stress plots 
corresponding to initial and optimal shapes are shown in Figure 5.28. A total 32.9% 
reduction in frontal-area of the fillet ABC has been achieved. 


5.2.5 A Ladle Hook Problem 


Application of the proposed shape optimization methodology for an elastic structural 
component is introduced as follows. 

Consider the design of a hook for lifting hot metal ladles with minimum weight of 150 
kN (hence the load on each hook is 75 kN). The dimensions of the critical elements of 
the hook (Figure 5.29 (a)) on the basis of stress, are presented by Dieter ([10]). These 
dimensions have been based on keeping the nominal stress at a level below 86.2 MPa. 
The remaining dimensions such as the bight section are set by engineering common sense. 




Figure 5.23 (a) Definition of a Fillet Problem with 

Displacement and Stress Boundary 
Conditions. 

(b) Boundary Element Mesh for a Fillet 
Problem. 
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Figure 5.26 (a) Cartesian Profiles for the Initial and 

Optimal Shapes of a Fillet. 

(b) Intrinsic Prifiles for the Initial and 
Optimal Shapes of a Fillet, 







nhapter 5. Case Studies Via Zero-order Method 


CARTESIAN SPACE 


• SHAPE DPTIMIZAnaN • 

Input Olrvctory 
« ir/utr/tKt/Jftkka/thcs UI 


mi«t Pobloa I 

(Via Z«ro’S3ird«u: UathodJ 


Uilng 

3-LIHCE 


Svlvct bilndcu using LEFT button, Prsss LEFT mou&t button for sscond cornsr 


Lflls cfnist 


JLfDs clntrin 
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of a Ladle Hook. 


(b) Intrinsic Profile for the Optimal Shape 


Of a Ladle Hook. 
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5,2.6 Sensitivity Analysis of An Elastic Ring problem 


The geometry belong to an optimal shape of an elastic ring has been considered here 
for sensitivity analysis. The effect of change in SDVs corresponding to an optimal shape 
has been studied via the sensitivity analysis procedure explained under Section 4.2. 

The change in frontal-area (i.e. merit function) and the area sensitivity with respect 
to variation in optimal shape design variables including; Ki and si. The variation in 
curvature Ki is denoted by Aki and arc length Si denoted by Asi. The change in 
frontal-area and sensitivities have been estimated for 20 different perturbation sizes. 

Figure 5.35 shows, the graph indicating frontal-area ($) versus changes in curvature 
(Aki). Figure 5.36 shows, the graph indicating area sensitivity versus change in curva- 
ture. Similarly Figures 5.37 and 5.38 shows, the variation in frontal-area ($) and area 
sensitivity with respect to the change in arc length (Asi). 
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Figure 5.36 Graph shows Area-Sensitivity Vs Change in 
Curvature ( A/ci )• 
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Chapter 6 


CASE STUDIES VIA 
FIRST-ORDER METHOD 


6.1 Software Development 


The computer programs which have been developed for non-linear equation solver, 
numerical integration, shape synthesis, automatic boundary element mesh generation 
and refinement, and Boundary Element Analysis for Design (BEAD) in Section 4.2 are 
also used in this method. In addition, the optimization module based on first-order 
method has been developed in FORTRAN 77 using the optimization algorithm outlined 
in Section 4.3. The multi- window grphics environment created as described in Chapter 4 
is also integrated with the proposed optimization cycle. In order to show the effectiveness 
of the proposed shape optimization methodology a numerical example of an elastic ring 
under diametral compression has successfully been implemented. 

6,2 Numerical Example 


The shape optimization methodology via First-order method described in Section 4.3 
is now applied to a mechanical component. The problem is characterized as plane-stress 
case and the intrinsic 3-LINCEs model {R-R-R model) is used in order to have shape 
changes. 

6.2.1 An Elastic Ring Problem 

This example shows an elastic ring under a single compressive load along its diameter, 
as showii in Figure 6.1(a). Taking the advantage of symmetry, one quarter of the ring 
is modelled, as shown in Figure 6.1, It is a case of plane-stress, and the objective is 
to minimize the weight (i.e. the minimization of frontal- area), while keeping the outer 
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boundaries defining the exterior shape of the ring unchanged. It is to be noted that the 
point B or C is allowed to move freely in the horizontal direction. A total of 100 nodes 
and 50 elements are used in the BEM analysis model. 

The formulation of the problem can be stated as follows. 

Minimize 
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and 


a;3 - 0.05625 < 0 


(Hi) 
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= 0 degrees 

and %j)z 

= 270 degrees 


( 6 . 8 ) 


Using the above formulation and shape optimization methodology via first-order 
method outlined in Section 4.3, an optimal shape of an elastic ring is achieved, cis 
shown in Figure 6.2. 

Figure 6.3 shows, the SUN-view graphics environment with cartesian and intrinsic 
geometries for the initial and the optimal shapes of the elastic ring. 




Figure 6.2 (a) Cartesian Profiles for the Initial and 

Optimal Shapes of An Elastic Ring (Via 

First-order Method). r i j 

(b) Intrinsic Profiles for the Initial and 

Optimal Shapes of An Elastic Ring (Via 
First-order Method). 
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Chapter 7 


CONCLUSIONS 


7.1 General Observations 


The goal of CAD approach for shape optimization is to implement an interactive 
redesign system that integrates optimization methods within a flexible and efficient 
computational tool, easy to use by design engineers. This interactive module is intended 
to create the missing link between BEM and CAD technologies and therefore it should 
constitute one of the key elements in the complex chain needed, to computerize the 
design cycle. 

In particular, the shape optimization methodology proposed in this research requires a 
multi-window graphics environment. This is due to the fact that, the intrinsic definition 
of the curve and/or surface does not give the real feel of the geometry. In order to see 
the immediate effect of change in intrinsic shape design variables in cartesian space, the 
multi-window concept has been developed, thus the handicap of the intrinsic geometry 
would be eliminated. The multi-window concept developed on SUN-workstation using 
SUN-view graphics and UNIX-platform play an important role in the CAD based shape 
optimization cycle proposed in this research. The same optimization cycle on other 
systems, with the same multi-window concept has been developed on HP 9000/800- 
series workstation with UNIX-platform using STARBASE graphics in C language. It 
has been found that the optimization cycle is faster on HP-workstation than SUN- 
workstation. It is due to the supercomputing capability of HP-workstation. Figure 7.1 
shows, such a graphics environment. 

During the past two decades, the display technologies and computational tools have 
undergone a rapid change. Moreover, a unique integration of computational and display 
devices has occured in the form of computer control display devices. As the means of 
computations and display have changed, it has become imperative to re-examine the 
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geometrical theories used for modeling design and manufacturing problems. It can be 
seen that research in the area of computational geometry during the past two decades 
has brought out many interesting techniques of geometric modeling such as, for example, 
solid modeling. The parametric form of modeling geometry has been found most attrac- 
tive from two view points. In the first place, it is suitable for rendering the geometry 
in a parametric form or computer-controlled display devices. Secondly, by changing the 
position of some control points or auxiliary points, it is possible to manipulate the shape 
of the geometrical entities such as a curve or a surface. Extensive use of the parametric 
form has also brought into focus some drawbacks especially from engineering design 
view point. Inparticular, when a geometrical shape is to be manipulated for better or 
improved performance, then the relationship between the performance indices and the 
locations of the control points is not explicit. It is in this context that the proposed 
work has explored the role of one other form of representation of geometry, namely, the 
intrinsic form. It has been found that if the shape of a geometrical entity is presented 
in an intrinsic form, then it can be used for display purposes as well as for engineering 
analysis, optimization and manufacturing considerations. The proposed methodology of 
shape optimization utilizes intrinsic definition of surface boundary by selecting a shape 
model in the intrinsic space. For the selected shape model, it is possible to characterize 
a set of shape design variables by assigning suitable numerical values to the shape de- 
sign variables. In this way, it is possible to define the geometry of a surface boundary 
completely. Any change in the values of design variables reflects in the corresponding 
change of the shape of the geometrical entity. 

The proposed work consists of developing robust and stable numerical schemes in order 
to carry out shape synthesis as well as shape optimization using concepts of intrinsic 
geometry and BEM analysis. The shape synthesis of 2-dimensional surface boundaries 
can be accomplished by solving the Equations (2.31) through (2.33) provided that the 
curvature is given as a function of the arc length. The curvature can be represented 
using one of the proposed shape models. The shape models consists of a number of linear 
curvature elements (n-LINCEs). As the number of LINCEs increases, the number of 
shape design variables also increases, and hence the complexity of the optimization 
problem is augmented. The LINCE models seem to work for engineering design and 
optimization problems. 

Any shape optimization problem requires a geometrical definition of shape. Tradition- 
ally, this has been accomplished by using parametric spline functions. The co-ordinates 
of the control vertices of the spline functions constitute the shape design variables. Such 
an approach invariably leads to a large number of design variables. In contrast the pro- 
posed investigation clearly indicates that the number of intrinsic shape design variables 
is definitely less than the number of shape design variables in the parametric approach. 

In addition, it has been concluded that the intrinsic geometry approach seems to be 
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able to provide direct control for manipulation of the shape of a curve. The next im- 
portant feature of intrinsic geometry shape definition is that, it provides an opportunity 
to include intrinsic properties like curvature and/or torsion in to the definition of the 
curve explicitly. These features of intrinsic geometry makes the shape synthesis module 
elegant and versatile, whereas the Bezier and B-spline representations do not provide 
the capability to control the intrinsic properties of a curve directly. The intrinsic ge- 
ometry approach, however needs a scheme of mapping the intrinsic description on to a 
cartesian frame. This requires solving a set of non-linear equations and numerical inte- 
gration schemes. In the present case Newton-Raphson modified method and Simpson’s 
Numerical integration scheme were used for numerical analysis work. At this stage, 
it is important to note that while selecting (2n — 2) free shape design variables it is 
necessary rather must out of three dependent design variables, the arc length (i.e. s,’s 
) variables can not be more than one. Otherwise, the numerical solution of non-linear 
equations (2.31) through (2.32) would not converge at all. This is due to the fact that 
the non-linear equations are too nonlinear in terms of arc length design variables (s,’s) 
than curvature design variables (k,’s). The robust numerical schemes adopted in shape 
synthesis eliminates totally the user’s intervention in the design cycle. 

The versatility of this methodology is that the integration between shape synthesis 
module and BEM analysis module has been coupled through an automatic BEM mesh 
generation and refinement procedure. This eliminates the designer’s intervention during 
the design process. Since BEM analysis is usually carried out in cartesian domain, the 
present investigation has shown that once the mapping of geometry from the intrinsic 
plane to the cartesian plane is completed, the stress analysis is relatively straightforward. 

It is however, to be noted, that the proposed separation of the shape synthesis model 
from the analysis model, allows the designer to make a judicious choice in defining the 
shape via simple forms and low number of design variables, while at the same time using 
a highly refined BEM analysis model for sufficient accuracy of end results. Therefore, 
the matter of overall efficiency of the performance of methodology, remains to a great 
extent in the hands of the user. 

The shape synthesis module and BEM analysis module has successfully integrated 
via zero-order (i.e. Exhaustive Search) and first-order (i.e. Interior penalty and hence 
DFP) Non-Linear programming (NLP) methods by taking proper precautions while for- 
mulating the objective function and other constraints of the problem. The proposed and 
very useful general CAD approach for shape optimal design of 2-D elastostatic prob- 
lems consists of defining the boundary shape by Linear Curvature Elements (LINCEs) 
model. It provides a flexible and compact geometric capability. The Boundary Element 
Method offers significant advantage over other engineering analysis tools like FEM, FDM 
etc. in shape optimization, especially for mesh data preparation. Using this method, 
the stress analysis could be rapidly performed, the regeneration and refinement of the 
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mesh elements representing the geometry was both straightforward and fast, and the 
optimization of the component shape could be carried out automatically. 

The stress sensitivity and area sensitivity computational procedures explained in Sec- 
tions 3.3 and 4.4 respectively, have sucessfully been implemented and the sensitivity 
plots are obtained. 

The overall proposed shape optimization methodology has been implemented success- 
fully for illustrative examples which include: a thick cylinder subjected to a constant 
internal pressure, an elastic ring under diametral compresion, a torque arm subjected to 
an axial and transverse loads, a fillet under an axial load and a ladle hook subjected to 
tensile loading via zero-order method and an elastic ring problem via first-order method. 
The presented method of shape optimization approach can be applied to a broad range 
of applications in order to show its effectiveness over the parametric approach. The 
intrinsic geometry based shape optimization enables a designer to fulfil not only the 
functional requirements of an optimal design but also satisfy the manufacturing re- 
quirements. The proposed methodology will definitely become an effective tool in any 
CAD/CAM environment. 

The following aspects must be considered while designing a computer software for the 
CAD approach to shape optimization. 

The first one deals with the interactive computer graphics. It is present in most of 
the CAD systems as well as in analysis programs. It is concerned with the geometric 
transformations in 2-D and 3-D space and their graphical display. In the present case 
mapping the geometrical models from intrinsic to cartesian domain and then viewing 
them in an interactive SUN-view graphical display. 

The second aspect is related to the geometric modeling. This poses a problem of 
conversion of models. The geometric models from CAD systems have to be converted 
in order to use in the analysis and optimization programs. For example, in the present 
research the shape has been defined in the intrinsic domain and the analysis and opti- 
mization are carried out in the cartesian domain. 

The third aspect is the problem identification and the use of the design variables from 
the early step of shape design to the final step of optimization. The design variables are 
expected to be the parameters in the shape design, analysis, and optimization programs. 

7.2 Future Work 


Although it has been tried to include the intrinsic shape design variables directly into 
the shape functions of 2-D elastostatic BEM analysis, owing to the time limitations of 
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the research it has become necessary to leave this attempt incomplete. As a future 
work it has got a clear scope to implement. The description of the future work in this 
direction include defining the shape functions in terms of intrinsic shape design variables 
explicitly rather than presenting them in terms of cartesian x~y co-ordinates. In addition 
the BEM analysis requires derivatives of the shape functions w.r.t the local co-ordinates 
of the boundary elements (i.e. w.r.t and the evaluation of Jacobian of transformation, 
d(0' Thus, the introduction of intrinsic geometric design variables directly into the 
shape function might pose numerical complexities and instabilities together with the 
programming difficulty. Robust numerical schemes must be selected in tackling this 
proposed approach of BEM analysis for shape optimization. The problem of imposing 
displacement and stress boundary conditions can be handled as usual case of BEM 
analysis. But at the same time the sensitivity computations become simpler. 

Eventhough the present work utilized curvature function [k;(s)] is assumed to be a 
set of Linear Curvature Elements (LINCEs), the work can be extended using a set of 
cubic spirals. A cubic spiral is a curve whose tangent direction is described by a cubic 
function of arc length (s). It can also be suggested that the same work can be extended 
using higher order curvature functions too. These approaches are expected to give more 
smooth profiles as compared to the one with Linear Curvature functions, at the cost 
of increased computational complexity. This proposed approach has got wide range of 
applications which include robot path planning and optimal design of trajectories etc. 
The other applications include designing of optimal cam shapes for jerk free cam-follower 
mechanisms. 

The work can also be extended in the direction of designing optimal surface patches 
using intrinsic curve definition and lofted or ruled surface concept. For example, opti- 
mal blended transition surface patch design between air-craft fuselage surface and wing 
surface, optimal surface patch design for joining square or rectangular cross-section with 
the circular cross-section of air-conditioning ducts etc. by imposing suitable geometric 
and strength constraints. 

Using the intrinsic definition of surface and the three-dimensional elastostatic BEM 
analysis, intrinsic geometry based shape optimization methodology can be developed in 
order to tackle 3-D problems. This widens the effective use of intrinsic geometry in the 
shape optimization process. 
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